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Abstract 

A toy model for thermal contact consists of a two-spin system, where each spin is flipped 
by a thermostat. The transition rates are determined from the modified detailed balance 
discussed in Ref.[l]. Generalized heat capacities, excess heats, the housekeeping entropy flow 
and the thermal conductivity are calculated. The joint probability distribution of the heat 
cumulated exchanges at any time is computed explicitly. We obtain the large deviation func- 
tion of heat transfer via a variety of approaches. In particular, by a saddle-point method 
performed accurately, we obtain the explicit expressions not only of the large deviation func- 
tion, but also of the amplitude prefactor, in the long-time probability density for the heat 
current. The following physical properties are discussed : the effects of typical time scales 
of the mesoscopic dynamics which do not appear in equilibrium statistical averages and the 
limit of purely energy dissipation towards a thermostat when its temperature goes to zero. 
We also derive some properties of the fluctuations in the two-spin system viewed as a thermal 
machine performing thermodynamical cycles. 
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1 Introduction 



1.1 Issues at stake 

Though a rather specialized topic within non-equilibrium physics, the problematic of thermal con- 
tact is in itself a vast subject, and one of foremost theoretical and practical interest. A theoretic, 
microscopic, understanding of non-equilibrium dynamics is still lacking today, even in the particu- 
lar case of heat exchanges. When compared with the most general non-equilibrium physics world, 
certain specific settings have some advantages and this is the case of thermal contact: when only 
quantities conserved by dynamics (energy for a thermal contact) are exchanged between baths 
and an out-of-equilibrium system, and different baths are in contact with different parts of this 
out-of-equilibrium system, the fact that the baths remain at equilibrium all along the experiment 
allows to keep a simple, thermodynamical interpretation of various physical quantities (usually 
entropy variations) related to the exchanges. 

There are at least two important microscopic versions of thermal contact one can have in mind. 

One is as an interface between two thermal baths, generically at different temperatures. In 
our real (three-dimensional) world, the most natural model geometry for the interface is a real or 
fictitious (two-dimensional) surface. Each side of the surface consists of atoms of one bath. Heat 
flows from the high temperature bath to the low temperature bath either via interactions among 
atoms sitting on each side of an immaterial interface (a case relevant when the baths consist of 
solid materials) , or via interactions of the atoms with a thin or structureless material interface (a 
case relevant for instance when the baths are gaseous, no matter is allowed to be exchanged, and 
the interface is a piston). One can of course generalize to more than two baths, in contact via an 
appropriate number of interfaces. 

Another version of thermal contact is as an extended piece of material, with specified physical 
properties, and such that two (or more) parts of its boundary are in contact with thermal baths. 
A typical example is a bar of metal whose extremities are maintained at different temperatures. 
Though the piece of material is not in thermodynamic equilibrium stricto sensu, it is often of 
practical value to assign a local temperature at each point of the bar, and the famous Fourier 
law states that, for an isotropic material sustaining small temperature gradients, the heat current 
is proportional to the temperature gradient at each point. This law is of course approximate 
and phenomenological. But finding a physically motivated model, even a crude one, for which 
the relation between local temperature and heat current can be derived from first principles, and 
compared with Fourier's law, is a major challenge in the field. This is true even when the piece of 
material is a homogeneous bar and only its longitudinal dimension has to be taken into account, 
its section being homogeneous with a good accuracy. 

Our aims in this article can be considered as very modest, especially when compared to the 
general issues raised by both versions of thermal contact alluded above. But before describing 
in detail our model, let us put it briefly in context and compare it to other, more or less similar 
approaches. 

A significant trend in out-of-equilibrium statistical mechanics has been the search for solvable 
models which could give some hints in the comprehension of out-of-equilibrium phenomena in the 
absence of any theoretical framework which would play the role of Gibbs' statistical ensemble 
theory for equilibrium states. Several kinds of models for heat conduction have been introduced. 
In some of them the two heat baths are connected by a system with deterministic dynamics, 
such as an anharmonic chain (Fermi-Pasta-Ulam model) (5) |3J |3J [5J |B] or a one-dimensional hard 
particle gas [TJ [5] . In other models the system which ensures the energy transfer from one energy 
reservoir to the other has a stochastic dynamics: it may be an Ising spin system (see for instance 
among others [9] [TO] HU [12] ) or a particle with Langevin stochastic dynamics [13], [14] . In the 
latter case the heat exchanges are described as the work performed by a force including a friction 
as well as a random noise component. The latter interpretation of heat has been proposed and 
investigated by Sekimoto |15| for ratchets models, and has been used again in the interpretation 
of the Hatano-Sasa identity [TB] as well as in the investigation of heat fluctuations in Brownian 
transducers [T7] . 
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In the first part of the present work [T|, referred to as paper I in the sequel, we have inves- 
tigated the generic statistical properties for experimentally measurable quantities in the case of 
thermal contact models with the following features: the system has a finite number of possible 
configurations, the heat exchanges are described as changes in the populations of energy levels, 
and the configurations evolve under a stochastic master equation with transition rates bound to 



obey the modified detailed balance restated in (2.3 1. We have shown how the latter relation arises 



from the existence of an underlying ergodic deterministic microscopic dynamics which conserves 
energy. 

In the present article, we concentrate on the crudest possible version of a thermal contact: 
two heat baths are connected via the smallest possible contact. So in the "thermal contact at an 
interface" image, we would replace a surface by a point. And in the "heat flow through a piece of 
material" image, the bar we consider has microscopic length, cutting any hope of understanding 
temperature gradients. Though we believe that the model has some interest with regard to these 
interpretations despite its disarming simplicity, we shall give in a moment a third view of thermal 
contact for which our model is maybe less caricatural. But before that, we make yet another 
simplifying assumption on the way heat is exchanged between the baths: heat bath 1 (resp. 2) 
can flip a dynamical variable <t\ (resp. 02) which can take only two distinct values. The energy 
£(01, cr 2 ) changes when a contact dynamical variable is flipped, and assuming an energy conserving 
dynamics, this means that some energy comes from, or is given to, the heat bath responsible for 
the flip. Without loss of generality, we may assume that the values taken by 01 and 02 belong 
to {—1, 1}, and we shall use the name spins for a\ and 02- In general, the energy £((71,172) for 
the contact dynamical variables could take 4 distinct values, but we even concentrate on the case 
£(01, (7 2 ) = ^(1 — o^o^Ae, where Ae > is the energy gap. In the language of spins, this means 
the absence of external magnetic fields. More abstractly, it implies a twofold symmetry. 

We shall concentrate on a description of the time evolution of interface states (01,02) by 
a Markov process, i.e. by a probabilistic description. But other approaches are possible. For 
instance, as explained in paper I , the motivation for our choices of transition rates comes from 
invoking an ergodicity argument for a deterministic discrete time evolution of the compound "heat 
bath 1 plus interface plus heat bath 2". 

As the toy model has only 4 states, solving it can be reduced in some sense to the diagonalization 
of a 4 x 4 matrix, and the twofold symmetry of the energy functional allows to reduce this task to 
the diagonalization of a pair of 2 x 2 matrices. However, this is not the end of the story, and this 
takes us to the third interpretation of the model. 

The third view of a thermal contact mentioned above is not microscopic but mesoscopic. We 
regard <j\ and a 2 as some relevant collective variables and £((71,172) as an effective energy. Then 
the system can be viewed, and analyzed, as a thermal machine. That is, our crude model keeps 
track of one (and maybe only one) interesting feature: the system can make cycles. Consider a 
sequence of flips in the interface, starting from the state (o-i,0 2 ): 

(01,02) -> (-0-1,0-2) 

t I 
(0-1,-0-2) <- (-0-1,-0-2) 

after which the interface has returned to its original state. Writing E\, E 2 for the initial energies 
in heat bath 1 and 2, the sequence translates into 

(E U E 2 ) -> {E x -at 0-2 Ae, E 2 ) 

I 

(Ei - 2o-i0 2 Ae, E 2 + oi0 2 Ae) <- (E\ - 0-i0 2 Ae, E 2 + oi0 2 Ae) 
i 

(Ei - 2<ri<j 2 Ae, E 2 + 20io- 2 Ae) 

i.e. an amount of heat 2ai02Ae has been transferred from heat bath 1 to heat bath 2. We use the 
term mesocopic (as opposed to macroscopic) for two reasons: first the dynamics at the interface 
is not deterministic, i.e. knowing (01,0-2) at some time does not allow to know its value in the 
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future and second (this is somehow a consequence though) there may be portions of time in which 
the net flow of heat is from the cold bath to the hot bath. 

When seen in this light, the model is already more interesting: the time evolution of the 
heat bath energies is a random walk in continuous tim^] a subject known to lead to a number 
of nontrivial mathematical problems, some of them having a direct physical relevance. And in- 
deed we shall concentrate mainly on the physics, with the aim of performing detailed analytical 
computations. 

Another interest of our specific solvable model is that it plays the role of a pedagogical example 
where the general statements are made very explicit. For instance, though the fluctuation relations 
entail a constraint upon large deviation functions, they do not allow to determine them. The 
analytical calculation of the large deviation functions may provide a deeper understanding in the 
information which they contain. We shall see that, within the model, the computation of large 
deviation functions for the energy variations in the baths can be remarkably simple or tricky, 
depending on the kind of techniques one uses. 

Finally, let us note that in the absence of any general framework for out-of-equilibrium statisti- 
cal mechanics, the formulae obtained for the solvable model can give a flavor of the physical effects. 
For instance, the time scales of the microscopic dynamics, which do not show off in equilibrium 
averages, play a role in out-of-equilibrium properties even at the macroscopic level. Moreover the 
model can be considered in the limit where the temperature of the cold thermostat vanishes ; then 
the strict dissipation of energy towards the zero-temperature gives rise to specific phenomena. 



1.2 Contents of the paper 

The results of the explicit analytical calculations for the solvable model where the system is reduced 
to two spins are the following. 

In the case where the spin system involves only two spins, the transition rate are determined by 



the modified detailed balance (2.3) up to the typical inverse times v a of spin flips by each thermal 
bath a. characterized by its temperature T a . For an Ising interaction between the two spins, the 
transition rate for the energy exchanges with one bath take a form similar to that introduced by 
Glauber |T5] in his investigation of the time-dependent statistics of the Ising model in contact 
with one thermal bath and which is enforced to evolve towards equilibrium in the infinite-time 
limit. Most of the time, our results will hold whatever the values of T\ and T2 are. However, it 
is sometimes convenient to know in which direction heat flows from one reservoir to the other, 
and then we shall always assume that T\ < T2. By symmetry, this induces no loss of generality 
anyway: the results for T± > T2 can be retrieved by permuting T± with T2 and v\ with v-i ■ 

The Non-Equilibrium Stationary State (NESS) of the model happens to have a very specific 
property : since the transition rates are invariant under the simultaneous flips of both spins, 
the configuration probability distribution in the NESS coincides with an equilibrium canonical 
distribution at some inverse temperature /?+. 

The linear and non-linear static responses are discussed (see section [3]). Generalized heat 
capacities are defined; their explicit expressions involve the typical inverse time scales i/ a 's of the 
heat exchange dynamics with the energy reservoirs a = 1,2. The v^s are also called kinetic 
parameters in the following. In the vicinity of equilibrium the mean heat current is proportional 
to the difference T2 — T\ between the bath temperatures ; then a linear thermal conductivity can 
be defined. When the system is far from equilibrium the mean heat current is a bounded function 
of the thermostat temperatures (saturation phenomenon) ; one can introduce a non-linear thermal 
conductivity which vanishes in the limit where the relative temperature difference goes to infinity. 
The housekeeping entropy flow and the excess mean heats are defined in terms of the measurable 
averages of the cumulated heats and they are explicitly calculated for the model. 

The joint probability distribution for received cumulated heats Q\ and Q2 is determined at any 
finite time and for any initial distribution probability through a generating function method (see 



1 And in two spatial dimensions. The sum of the two coordinates can only take a finite number of values, but 
the waiting times and some correlations prevent from concentrating only on one component. 
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section [4]). Other distribution probabilities are then derived from its expression (4.42), and the 
explicit results are summarized in subsection |4.2| The results are given in terms of two integrals in 
the complex plane. The system not only obeys the finite-time symmetry enforced by the modified 
detailed balance for the ratio of the probabilities to measure some given heat amounts Q± and Q2 
or their opposite values when the system is initially prepared in an equilibrium state, but it also 
satisfies another finite-time symmetry specific to the model for the ratio of the same probabilities 
when the system has any initial distribution probability. The latter fluctuation relation involves 
the initial probability for the sign of the spins product, probably because, by virtue of energy 
conservation, the values of this sign in the final and initial states are related to the sum Q2 + Qi 
by the constraint (4.51. 

The cumulants for the cumulated heat Q2 are studied in section [5] from the characteristic 
function of the probability density for 0,2- The relation between the characteristic function of a 
probability density II(<2; t) with the generating function for the probability function P(Q; t) when 
the variable Q can take only discrete values is recalled in subsubsection [5~1~T] The explicit formulae 
for the first four cumulants per unit time in the infinite-time limit are given in (5.151. Even at 
equilibrium the cumulants arc not those of a Gaussian. 

The large deviation function for the cumulated heat current Q2 jt is calculated by three different 
methods (section [6]) : from the Gartner-Ellis theorem (subsection 6.1 ), from a saddle-point method 
(subsection |6.2[) and from Laplace's method on a discrete sum (subsection 6.3). The second and 



third methods rely explicitly on the discrete nature of heat exchanges in the model and on an ad 
hoc definition of large deviation functions, but they allow to compute subdominant contributions 
as well. The first method is straightforward, one just has to check that the general applicability 
hypotheses (recalled in detail below) are fulfilled and this is easy in our case. The third method 
is also simple because it deals with a sum of nonnegative terms, so no compensation is possible. 
The saddle point method however is remarkably tricky in our case, for reasons that we shall detail 
below. The expressions in terms of various parameter sets are given in (6.10 1 and (6.16). In order 
to readily obtain the large-deviation function in the case where the temperature of the colder 
bath vanishes, its expressions for positive and negative currents are explicitly distinguished in 
(|6TT3|)-([6TT4]). 

The limit where the kinetic parameter of one thermostat becomes infinitely large with respect 
to the kinetic parameter of the other thermostat is studied in section [7] In this limit the stationary 
distribution of the spins is the equilibrium canonical probability at the temperature of the "fast" 
heat bath while the typical inverse time scale in the mean instantaneous heat current is the kinetic 
parameter ^ s i ow of the "slow" heat bath. The probability distribution for the heat amount received 
from the slow thermostat, Q s iow, at any finite time t is that of an asymmetric random walk with 
the inverse time scale v s \ ow . As a consequence the probability distribution for Q s iow obeys a 
fluctuation relation at any finite time (see (7.20 1). The probability distributions of <7i<72 and Q s i ow 
are independent from each other, and this mean-field property is interpreted as a kinetic effect in 
the considered limit. The very simple forms of the long-time cumulants per unit time are given. 
The long-time distribution of the cumulated heat current is exhibited : it vanishes exponentially 
fast over a time-scale given by the inverse of the large-deviation function (7.33 1 with an amplitude 
which is explicitly calculated. 

In the limit where the temperature of the colder thermostat vanishes (section [8| the microre- 
versibility is broken but the system still reaches a stationary state where all configurations have 
a non- vanishing weight, because the Markov matrix is still irreducible. The large deviation func- 
tion is expressed in (8.9 1. In the limit where ^2/^1 +00 (resp. v\jv<i — > +00) the probability 
distribution for the heat amount —Q\ (resp. Q2) at any finite time t becomes a Poisson process, 
because the zero-temperature thermostat can only absorb energy (pure dissipation towards the 
zero-temperature bath) . Again the very simple forms of the infinite-time cumulants per unit time 
are given, as well as the large deviation function (8.21 1. 

The last section is devoted to a probabilistic study of the system seen as a mesoscopic thermal 
machine, with an average heat flow from the hot reservoir to the cold reservoir, but also fluctuations 
around the average which we try to quantify. We compute the probability for the thermal machine 
to work backwards, and the law of the fluctuations of the time it takes to the machine to do one 
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cycle. We present the argument for a system slightly more general than the two-spin system, 
because the computations and their meaning are more transparent this way, and then apply the 
formulae to the two-spin system. 



2 Model 

The physical system we deal with in this article is a toy model of thermal contact, consisting of two 
heat baths, generically at different temperatures, put indirectly in contact via a small subsystem 
made of two interacting Ising spins u\ and o 2 . Each spin a a , a = 1, 2 is in contact with a single 
bath denoted by a. We aim at a statistical description, where the details of what happens in 
the heat baths is not observed, but only the evolution of the two spins, i.e. of the configuration 
C = (oi, <r 2 ). We assume that this evolution is described by a Markov process (in continuous time) 
with transition rate (C'|W|C) from configuration C to configuration C . 

As the system is out of equilibrium, the form of W is not a direct consequence of known physical 
laws, and it is unclear whether a nature-given preferred choice exists. So we start with a purely 
technical and down to earth description of our choice for the transition rates, that we shall use for 
all later explicit computations. The general principles and steps that guided us to the modified 
detailed balance that the transition rates must obey have been given in paper I . The main ideas 
are the following. 



2.1 Constraints upon transition rates arising from microscopic discrete 
ergodic energy-conserving dynamics 

As usual, we view a heat bath as an ideal limit of some large but finite system. So the system 
we describe is obtained via a limiting procedure from a large system made of two large parts and 
a small one, which is reduced to the two Ising spins oi and o 2 , each one directly in contact with 
one of the large parts. 

We expect that in this limit many details become irrelevant, so we assume for the sake of the 
argument that the degrees of freedom in the large parts are discrete. 

As in classical statistical mechanics, we take the viewpoint that the statistical description 
of 01,02 is an effective mesoscopic description arising from a deterministic, energy-conserving 
dynamics for the whole system. With discrete variables, there is no general definition of time 
reversal invariance, but we impose that the dynamics is ergodic. 

We also want the dynamics to reflect the fact that the two large parts interact only indirectly: 
there is an interaction energy £(01, 02) between the two spins and the spin o a is flipped thanks to 
energy exchanges with the large part a (a = 1, 2). Defining the operator ¥ a as the operator flipping 
the spin o a while leaving the other spin unchanged (e.g Fi(oi,o 2 ) = (—01,02)), in this process 
the energy of the large part a is changed from E a to E' a according to the energy conservation law 

K _ Ea= {-[£(C')-£(C)\ if C' = ¥ a C 
10 otherwise, 

while the energy of the other large part is unchanged. 

As shown in paper I , when the large parts are described at a statistical level and in a tran- 
sient regime where the large parts are described in the thermodynamic limit, the transition rate 
(C'|W|C) from configuration C to configuration C obeys three constraints: first the graph associ- 
ated with the transition rates is connected; second there is microscopic reversibility for any couple 
of configurations (C,C), 

(C'\y\C)^0 (C|W|C')^0; (2.2) 

third the ratio of transition rates obeys the so-called modified detailed balance (MDB), 

for C'-FC ( C 'I W I C ) - e -Pa[s(C)-S(C)] (o o\ 
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We remind the reader that the latter relation is also referred to in the litterature as the "generalized 
detailed balance". 

2.2 Determination of transition rates 

The transition rates are non-zero only if the initial and final configurations C and C differ only 
by the flip of one spin : (C'|W n |C) = unless C = F a C. Since a a can take only the two values +1 
and — 1, the transition rate where a a is flipped takes the generic form 

(F a C|W|C) = ^^[l- ( r Q r Q ( ( 7 b )]. (2.4) 

The four parameters v ai T a , a — 1, 2, are a priori arbitrary except that the z/s are > and the F's 
are of absolute value < 1. 

Taking for simplicity an interaction energy between the spins 

£(a 1 ,a 2 )= 1 ~^ ia2 Ae, (2.5) 
where Ae > is the energy gap between the two energy levels, one gets from the modified detailed 



balance in the form (2.3) that 



a condition similar to the one obtained by Glauber [18] in the equilibrium case. As e 2x = Manila 



the generic form (2.4 1 of (F a C|W|C) has to satisfy 



(¥ a CMC) = V -^- [1 - o x o 2la \ (2.7) 

with 

7o = tanh LSb^V (2.8) 



If Pi and /?2 are finite, < 71 < 1 and < 72 < 1, and the microscopic reversibility condition \2.2\ 
is also satisfied. Without loss of generality we could, and will sometimes, assume that T\ < T 2 . 
Then 71 > 72. 

For the sake of simplicity, in the following we assume that v a depends only on the properties 
of the thermostat and not on the value of a\,- (This choice enforces the equality between the 
transition rate from ((71,172) and that from (— ctj.,— 02), which are two configurations with the 
same energy.) Apart from simplicity, we have no convincing argument that this should be THE 
nature-given preferred choice. Anyway, we write 

(F„C|W|C) = ^[l-(7l(727a]- (2.9) 

This ends the argument explaining our choice of transition rates and gives a physical interpretation 
of the parameters: 7 a is formed with the energy scale in the two-spin system and the temperature 
of the bath a, while v a describes a rate at which the bath a attempts to flip spin a a . 

We notice that, though the transition rate expressions have been derived from hypotheses 



implying the microscopic reversibility (2.2), these expressions still make sense if P2 < Pi — +00. 



(The limit /3± — > +00 where microscopic reversibility is broken is discussed in section |8j) 
Moreover, even if fi\ — +00, the Markov matrix M defined by 

(CW) = I™ C) , ^ C '/ C (2.10) 

is irreducible, namely any configuration C can be reached by a succession of jumps with non-zero 
transition rates from any configuration C. 
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3 Non Equilibrium Stationary State (NESS) as a canonical 
distribution with an effective temperature 

3.1 Stationary state distribution 

The master equation which rules the evolution of the probability P(C;t) can be written in terms 
of the Markov matrix M defined in (2.101 as 



dP(C;t) 
dt 



(3.1) 



In the basis where the probability P(ai,(T2;t) is represented by the column vector 



\P(t)) 



the matrix M takes the form 



V\ + h> 2 



-1+7* 


^2(1 - 72) 

\l7i(l - 71) 



fp(+ 


+;t)\ 




P(- 


-,t) 




P(+ 


-,t) 




\P(- 


+;t)J 







^2(1 + 72) 


+ 7* 


^1(1 + 71) 


-7i) 




1-7* 


-72) 








(3.2) 



^2(1 + 72' 


-1-7* 



(3.3) 



In the latter equation we have introduced the dimensionless inverse time scales 

V„ = — — — for a = {1, 2}, 



v\ + v 2 



and we have set 



7* 



f l7l + ^272- 



(3.4) 



(3.5) 



The Markov matrix ¥\ is irreducible (even if 7 = 1, namely T± — > 0): for any pair of configura- 
tions C and C, there exists a succession of spin flips, with non-zero transition rates, which allows 
to make the system evolve from C to C' . Henceforth, according to the Perron- Frobenius theorem 
there exists a single stationary state distribution P s t(C) and it is nonzero for every configuration 
C. 

Moreover, since the system is made of two discrete variables which can take only the values ±1 
and since the transition rates are invariant under the simultaneous flips o\ — > —a± and 02 — &2, 
the stationary distribution P s t(ci, 02) takes the form P st (<7i, 02) = a + da\<T2- Indeed, the generic 
form of P(<7i,o"2) reads P(a\,/J2) — a + bo\ + ca 2 + do\<Ji. On the other hand, the invariance 
of the transition rates under the simultaneous flips &\ — > —<J\ and 02 — > — C2 entails that if 
a + b<7\ + C02 + do\Oi is a stationary solution, a — bo\ — co-i + da±a 2 is also a stationary solution. 
But, since M is irreducible, there is only one stationary solution, so that b = c = 0. By solving 



explicitly the master equation (3.1) and using the normalization of a probability distribution, the 
stationary solution proves to be 



P s t (0-1,0-2) 



1 



7*°-io-2j 



(3.6) 



where 7* is defined in (3.5 1 



The stationary distribution of the model has the following remarkable property: it coincides 
with some equilibrium distribution. More precisely, the stationary state distribution is equal to 
the canonical state distribution at the effective inverse temperature /3* 



Pst (o-i,o- 2 ) = (0-1,0-2), 



(3.7) 
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where /?* is determined by the relation 



7* = tanh(&~J, (3.8) 

and 

e -/8£(C) 

^can(^) = ~~ Z~(j3)~ ' (3-9) 

where Z(/3) is the canonical partition function at the inverse temperature /3, Z{fi) — ^ c e~^ £ ^ c \ 
We notice that the canonical form for the state distribution implies that /3* obeys the canonical 
ensemble relation which is equivalent to the definition of the inverse temperature in the micro- 
canonical ensemble, namely 

e. = d ~™. (-») 

where (£) s t = J2c £(C)Pst(C) is the stationary mean value of the energy and S SG [P s t] is the value 
of the dimensionless Shannon-Gibbs entropy in the stationary state. The dimensionless Shannon- 
Gibbs entropy (where the Boltzmann constant is set equal to 1) is defined from the configuration 
probability distribution P(C; t) as 

S SG [P(t)] = -J2 P(C; t) lnP(C; t). (3.11) 

c 

Its evolution has been discussed in paper I . 



3.2 Linear static response to a variation of some external parameter 

In the present section we consider the static linear response of some observable O to a change of 
some external parameter, namely the inverse temperature f3 a or the typical inverse time scale v a 
of the bath a, with a = 1, 2. 

In the protocols for the study of static linear response, the system is prepared in some stationary 
state at time to = CP and the external parameters are instantaneously changed by infinitesimal 
amounts at time t = 0. Then, in the infinite time limit, the system reaches another stationary 
state corresponding to the new values of the external parameters. 



3.2.1 Relation with static correlations for a "canonical" NESS 



Since the nonequilibrium stationary distribution given by (3.6 1 involves only one parameter, 
namely /3+, the linear response coefficient d(0) s t/dg cx t for the mean value of an observable O 
in the stationary distribution when some external parameter g cxt is varied is proportional to 



d(0) st /d^, namely d(0} st /dg ext = (d^/dg ext ) x (d(0} st /d/3+). Moreover, by virtue of the 
stationary distribution is the canonical distribution at the inverse temperature /3*. Henceforth the 
coefficient d(0) s t/df3+ is merely opposite to the correlation between O and the energy £ according 
to the canonical equilibrium identity 

= - [(OS)tn - (0)tn(£)tn] ■ (3-12) 

where (0)^* D denotes an average with respect to the canonical distribution P^* n . As a result, 
the relation valid for responses to the variation of any external parameter in the nonequilibrium 
stationary state reads 

= [{OSU - (0) st (£) st ] . (3.13) 



9 



3.2.2 Dependance of the mean energy upon the time scales of the microscopic dy- 
namics 



The main difference between the response of the mean energy in non-equilibrium and equilibrium 
states arises for the response to a variation of the time scales of the microscopic dynamics which 
rules the heat exchanges with the baths. When j3\ = (3 2 the equilibrium mean energy (£) cq = {£)cin 
depends only on the thermodynamic temperature common to both baths. On the contrary, in the 
non-equilibrium case the stationary mean energy {£) s t does also depend on both inverse time 
scales v\ and v 2 . Indeed, since the stationary probability corresponds to the effective canonical 



distribution (3.7 1, the stationary mean energy reads 



(£) st = (1 - 7*) ^ = (1 - v xlx - V 2l2 ) ^. (3.14) 

Changing v a means changing the physical connection between the thermal bath a and the spin 
system. The linear response of the stationary energy associated with a variation of the inverse 
time scale v a is determined by the coefficient 

~di^~ = ~Jy^+v^ la ^2 fo r{a,6} = {l,2}. (3.15) 

3.2.3 Stationary mean energy and generalized heat capacities 

The heat capacity C oq is a measurable quantity defined as the ratio 

C cq (T) = ^, (3.16) 

where (6Q) is the mean heat amount received by the system in transformations which involve only 
heat transfers and make the system go from an equilibrium state at temperature T to another 
equilibrium state at temperature T + dT, while all other thermodynamic parameters which deter- 
mine the equilibrium state are kept constant. ( (5Q) = lim^ +00 (Q) t in the protocol mentioned 
in the introduction of the section.) According to the energy conservation, (SQ) = {£)l+ dT - (£ )J q 
and the heat capacity is related to a partial derivative of the equilibrium mean energy, 



r;, ( rn = ^ = -/3 2 ^. cud 



When the system is in a stationary non-equilibrium state induced by thermal contact with two 
heat reservoirs at respective temperatures T\ and T 2 , we can introduce measurable heat capacities 
by similar definitions. When the temperature T\ of the thermal bath 1 is changed by dT\, while the 
temperature T 2 of the thermal bath 2 is kept fixed, and when the system evolves from a stationary 
state to another one only by heat transfers, then the generalized heat capacity is defined as 

cl t 11 (T 1 ,T 2) = ^ ( ^+ g2)) . (3.18) 

According to conservation energy, S(Qi + Q 2 ) = (£) s t 1: 2 — (£) s i' 2 an d the heat capacity is 
related to a partial derivative of the stationary mean energy 



cll\T 1 ,T 2 )= d{£) ^ T2 



(3.19) 



T-2 



In the present model the expression ( 3.14[ ) of the stationary mean energy takes the very specific 
form 

(£) st = u 1 (£)^+V 2 (£)^. (3.20) 



10 



Indeed the relation v\ + v 2 = 1 and the expression of the equilibrium mean energy at the inverse 
temperature /?, 



Ac 



when /3i = /3 2 = P, 



(3.21) 



allow to rewrite the mean energy expression (3.14) in the non-equilibrium stationary state in the 
form ( |3.20[ ). By virtue of the specific decomposition (3.20 1 of the mean energy, the heat capacities 
C% ] (Ti,T 2 )'a read 

C^ ] (T 1 ,T 2 )=V a C cq (T a ) with a = {1,2}, (3.22) 



where, according to the relation (3.171 and the expression (3.211 of (£) 



C e q(T a ) 



1 - tanh 



(3 a Ae 



(3.23) 



More generally, when the temperatures T± and T 2 of both thermostats are varied independently 

C 8 [ t 1] m , T 2 )dT x + CH ] (Tj , T 2 )dT 2 . (3.24) 



If Ti and T 2 are increased by the same infinitesimal quantity dT the corresponding heat capacity, 
defined as C st (T u T 2 ) = 5{Qi + Q 2 )/dT is equal to the sum C$(Ti,T 2 ) + Cf ] (Ti,T 2 ). For the 
present model C st (7i,T 2 ) = ViC^T^ +V 2 C cq (T 2 ). In the limit where T X =T 2 = T, by virtue of 
the relation Vi + V 2 = 1, we retrieve the equilibrium heat capacity C eq (T), as it should be. 

3.2.4 Stationary heat current and linear thermal conductivity 

The instantaneous heat current j a (C) received from the heat bath a when the system jumps out 
of the configuration C has been defined in paper I as 



j a (C) = JSqa (C) = J2( c 'M c ) s Qa(C 4- C), 



(3.25) 



where 8q a {C C) is the heat received from thermal bath a when the system evolves from config- 
uration C to configuration C = ¥ a C, where F a is the flip caused by thermal bath a, namely 



Sq a (C <- C) = [£(C) - £{C)\ if C = ¥ a C 
Sq a (C C) = otherwise. 



(3.26) 



In the present model j 2 (ai, a 2 ) = [£(ai,-a 2 ) —£(ai,a 2 )} (ax, —ct 2 \W\o-l, a 2 ) = v 2 \o\o 2 -72] (Ae/2). 
In the stationary state the mean energy is constant so that the mean currents received from 



both baths cancel, (ji) s t + (j2)st = 0. For the stationary state probability distribution (3.6), one 
has (aia 2 ) st 



7* 



f 171 + ^272 and 



(j2)st = viv 2 (71 - 72) 
where 71 — 72 may be rewritten as 



(vi + i> 2 )Ae 



71 — 72 = tanh 



l-tanh[M^ t anh^ Ae 



(3.27) 



(3.28) 



When T\ < T 2 , (j 2 ) s t > 0, as it should : the mean heat current flows from the hot bath to the 
cold bath. Note that {j 2 )st is a bounded function of T\ and T 2 . Thus, in the generic case (j2) s t 
is not proportional to the bath temperatures difference T 2 —T\. As for any system, the linear 
dependence upon T 2 — T\ (or /3i — f3 2 ) appears in the limit where (/?i — (3 2 )Ae -c 1. In the high 
temperature regime where both /3iAe <C 1 and (3 2 Ae <C 1, the condition (j3i — {3 2 )Ae < 1 is 
satisfied and {j 2 ) s t is proportional to T 2 — T\. 
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When fa = fa the system is at equilibrium and (j 2 ) C q — 0. Moreover, as shown in paper I, 
the partial derivatives of the current obey the generic symmetry 



d{j2) s t 



dfa 



2/st 



(fa = fafa=fa = ~ 



(fa=fafa=fa. 



(3.29) 



This property can also be checked from the expression (3.27 1 of (j 2 ). It entails that, d(j 2 ) s t / dT% \t 2 {T, T) = 
-d(h) at /dT 2 \ Tl (T,T) and subsequently 



IJ \ (rp rp N <9(j2)st 

w2)st , ~ (T 2 - Ti) 

(T U T 2 )^(T,T) 8T 2 



(T,T). 



(3.30) 



In other words, when Ti and T2 independently tend to the same value T, at first order in the 
independent variables T\ — T and T 2 — T the ratio (j 2 ) s t/{T 2 — T\) depends on T but is independent 
of the ways Ti — T and T 2 — T vanish. 

As a consequence, for a non-equilibrium stationary state near equilibrium, namely when the 
temperature difference between the thermostats is such that (fa — fa) Ae <§; 1, one can define the 
thermal conductivity as 



lim 



(h) s 



(Ti,T 2 )-KT,T) T 2 - Tx 



dT 2 



(T,T). 



'A 



From (3.271 we get the expression for the thermal conductivity, 

( fa\e 



^1^2 
V\ + v 2 



1 - tanh z 



V 2 



(3.31) 



(3.32) 



We remind the reader that the thermal conductivity, which is a positive transport coefficient, is 
related to the kinetic coefficient (also called Onsager coefficient) introduced in phenomenological 
irreversible thermodynamics as 

L= lim (3.33) 

where, as discussed in paper I, the thermodynamic force !F can be defined from the station- 
ary entropy production rate, which is opposite to the exchange entropy flow, di nt S SG /dt\ st = 
—dexchS/dt\ st , through the relation 



dt 



*^exch $ 

dt 



^(32 



(3.34) 



when there is only one independent mean instantaneous current. As shown in paper I , in the case 
of the thermal contact F = fa — fa. Therefore the relation between the kinetic coefficient and the 
thermal conductivity defined in (3.311 reads 



L 



«th 
fa ' 



(3.35) 



Now we compare the results about the linear static response in non-equilibrium stationary 
states which are either in the vicinity of equilibrium or far away from equilibrium. When the 
system is far from equilibrium, namely when (fa — fa) Ae ^ 1, (3.271 leads to 



d(32) s 



dfa 

9{j 2 ) a t 
dfa 



Pi 



V\V2 
V\ + V 2 



V\V 2 

V\ + v 2 



1 — tanh 



1 — tanh 



/3 2 Ae 



faAe 



Ae 
~2~ 



Ae 
~2~ 



(3.36) 
(3.37) 
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The linear response coefficients <9(j 2 ) s t/<9T 2 | Ti and d(j2)gt/dT\\ T are no more opposite to each 
other. As a consequence, when T\ and T 2 are varied independently, the corresponding variation 
of the stationary mean instantaneous current {j'i)ft a t first order reads 



<*>]?' 



\j2) s t 



(T[,T^(T U T 2 ) 



Irril rrt \ ^W^/st 1 (rnl rri \ w \J2/ s t 

T 2 



. (3.38) 



'A 



The latter variation depends not only on Ti, T 2 and the variation of the temperature difference 
(T{ — T 2 ) — (Ti — T 2 ) but also on the way in which T[ and are varied around the given values 
T\ and T 2 . 



3.3 Non-linear static response in the NESS 
3.3.1 Non-linear thermal conductivity 



MA) 



When the system is far from equilibrium, instead of introducing the linear response (j2/ s t 
O'a)!* 1 '^ 2 ' w ith (/3[— /?i)Ae -C 1 and (/3 2 — /?2)Ae -C 1 (and the associated linear response coefficients 



9(i2)st 1 '^ 2 V^a)) one ma y rather consider a non-linear thermal conductivity defined as 



„nlir 



(J2>i 



T2-T1 



From (3.271 



„nlin 



1 



tanh 



fe)Ae 



TiT 2 ^1 + iaj 



(Pi 



-fe)Ae 
2 



1 — tanh 



tanh 



/3 2 Ae 



Ae 
~2~ 



(3.39) 



(3.40) 



According to the expression (3.27 1, (j 2 ) s t is a bounded function of Jf and T 2 , so that k"^ 111 vanishes 
when T 2 — Ti becomes very large with respect to either T\ or T 2 . 

We also notice that when both thermostats are at very high temperature, namely when fii Ae <C 
1 and /3 2 Ae <C 1, (j 2 ) s t is proportional to /3%— /3 2 with a coefficient independent of the temperatures. 
As a consequence, the partial derivatives d(j 2 )gf 1 '^ 2 '/<9/?i and c*^)^ 1 '^ 2 ' /<9/3 2 are opposite to each 



other, as in the symmetry property (3.291 in the very vicinity of the equilibrium limit. Then the 
difference (3.381 is proportional to the difference (f3[ — j3i) — (f3' 2 



/3iAe->0,/3 2 Ae-!-0 



ft), 

Pi) ~ " ft) 



90' 



2/ St 



dft 



(3.41) 



Besides, the thermal conductivity (3.40) behaves as 

1 



/3iAe-)-0 TiT 2 V\ 
/3 2 Ae->0 



^2 



Ae 
~2~ 



(3.42) 



3.3.2 Housekeeping entropy flow and mean excess heats 

In the long-time limit, whatever the initial configuration probability Pq may be, the system 
reaches a stationary state where the Markovian stochastic dynamics enforces that the cumu- 
lated heats received from each thermostat, namely the random variables Qi(t) and Q 2 (t), have 
averages (Qi{t))p and (Q,2(t))p which both grow linearly in time with opposite coefficients, 



-<Ci(*)>fb 



t— >+oo 



'Mt)h 



t—y+oo 



Then 



lim 

i— >+oo 



Pl{Ql{t)) Po +fi2(Q2{t)) i 

t 



-(A 



ft)(j 2 ) st 



(3.43) 
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where the stationary exchange entropy flow appears by virtue of (3.341. Meanwhile the sum 



Qi{t) + Q 2 (t) remains bounded at any time and its average tends to the heat amount corresponding 
to the mean energy difference between the final and initial stationary states, 



lim (Q 1 (t) + Q 2 (t)) Po ^(£) [ s ^ M ^(£) Po . 



(3.44) 



In the phenomenological framework of steady state thermodynamics [TU] , when work is supplied 
to the system, the total heat given to the system is usually expressed as the sum of an "excess" 
heat Q cxc associated with the energy exchange during transitions between two different steady 
states and a "housekeeping" heat Qhk associated with the energy supplied to maintain the system 
in the NESS reached in the long-time limit. These two heat amounts have been discussed for a 
system in contact with only one thermal bath and submitted to a time-dependent external force 
which is described by Langevin dynamics |16| |2"U], 

By analogy, with the standard sign convention, we may introduce a "housekeeping" entropy 
flow supplied to the system which can be measured as the asymptotic behavior 



Chk[-P s t] = - lim 



fa{Qi{t)) Po +fa{Q 2 {t)) Po 



(3.45) 



and which, by virtue of ( 3.43| coincides with the opposite of the stationary exchange entropy flow, 
namely with the stationary entropy production rate (see ( 3.34 1 ) 



0hk[-P s t] = - 



^exch 

dt 



S 



dint S 

dt 



(3.46) 



From the explicit expression (3.27) of the mean instantaneaous heat current we obtain the expres- 



sion for the housekeeping entropy flow (3.45) 



fhk[-Pst] = 



V\V2 
V\ + v 2 



in - 72) {fa - fa) 



Ac 



(3.47) 



When the system is prepared in a stationary state by thermal contact with heat reservoirs 
at the inverse temperatures /3° and fa] respectively, then (£}p = (£)st an d the difference in 



fl3~44) becomes equal to (f)^ 1 '^ 2 



With the standard convention, we introduce the 



"excess" heats given to the system Q c 



|[/3l,/3: 



„oi with a = 1, 2 which can be measured as 

1 •P 2 i 



Q 



I [01.02] 

1 •■'■<<*\[/3°,p°] 



lim 

t— >-\-oo 



(Q»(*)>, 



,LB?,/3§ 



Then, by virtue of the stationary condition (ji) s t = — (j 2 ) s ti the equality (3.441 becomes 



n lift" 82 ] 

Vcxc.H^o^O] 



^oxc,2|^O j/3 0] 



(3.48) 



(3.49) 



In the linear response regime where the relative differences (T\ — T®)/T® and (T 2 — T 2 )jT® are 



infinitesimal, by virtue of the definition (3.191 of the generalized heat capacities (Ti,T 2 ), with 
o = l,2, 

n i[T 1 +dTi,T 2 +dT 2 ] n I pi +dTi, T 2 +dT 2 ] /rp r UT W 2 1 ("T T 1 \JT ft rnl 

— Wexc,l\[ TliT2 ) - Voxc,2|[ Tl T ,i -? O gt (±x, l 2 )al\ + O st (li, l 2 )al 2 . (d.OUJ 



lpi,T; 



The excess heats Qexc olfao'floi s defined in (|3.48| are explicitly calculated in subsection 



4.3 



Joint probability distribution for heat cumulated exchanges 
at finite time in the model 



Instead of studying the evolution of the probability distribution P(C\t) of the spins config- 
uration C — {oi,<j 2 ), we address directly the evolution of the joint probability distribution 
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P (C'\Qi, 0,2, t\C) for the cumulated heats Q\ and Q2 received from the thermal baths 1 and 
2 during a time t when the system is in configuration C = (a%, <j 2 ) at time t Q = and in configu- 
ration C — {cr'ijCr^) a t time t. In order to obtain results which hold as generally as possible, the 
initial probability distribution for configurations is not assumed to have the same symmetry under 
simultaneous spin flips as the stationary distribution. 

Since the two-spin system has only two energy levels separated by the energy gap Ae, the 
cumulated heats Q a are integer multiples of Ae and we set 

Qi = -niAe and Q 2 = n 2 Ae. (4.1) 

The minus sign in the definition of Qi is introduced for the sake of conveniency, because the mean 
instantaneous heat currents (ji) s t and (j2)st in the stationary state are opposite to each other. In 
other words, niAe is the amount of heat dissipated towards the heat bath 1, while r^Ae is the 
amount of heat received from the heat bath 2. With these notations P(Cf\Qi, Q 2 ,t\Co) can be 
written as a matrix element of some evolution operator U(m, n 2 ; t) 

P(C'\Qi,Q 2 ,t\C) = (a' 1 ,a 2 \V(n 1 ,n 2 ;t)\a 1 ,a 2 ). (4.2) 



4.1 Explicit calculations 

4.1.1 Constraint from energy conservation 



According to the expression (2.5) for the interaction energy between the two spins, the energy 
difference between the final and the initial configurations reads 

S(a' 1 ,a' 2 )-S(a 1 ,a 2 ) = ai(72 - aW2 Ae, (4.3) 



and it can take only three values 0, +Ae and — Ae. On the other hand, according to (4.1|, 
Qi + Q2 = {n 2 - ni)Ae, namely 

Qi + 0,2 — An x Ae where An = n 2 —n\. (4.4) 

The energy conservation entails that the energy variation of the two-spin system is equal to the 
sum of the heat amounts received from the thermostats: £(a[,a' 2 ) — £(ai,cr 2 ) = Qi + Q 2 - As 
a consequence the correspondence between the total amount of received heat and the couple of 
initial and final states reads 

An = = (T1CT2 (4.5) 

(An) 2 = 1 ■<=> o\a 2 = An and crjcr 2 = —An. 

Therefore it is convenient to introduce the decomposition 

U(ni,n 2 ;t)= <5» l2: n 1 +An l U(ni, "i + An; t). (4.6) 

An=0, + 1,-1 



In the basis {(+,+),(—,—),(+,—),( — h)} already used in ( 3. 2 1 the correspondence (4.5 1 enforces 
that U (ni ,n 2 ',t) can be decomposed into three 4x4 matrices 

U = ^[An=0] + U[A«=+1] + U [An =-l] ( 4 -7) 

with 

( k 0\ nl /O 0\ ,, f<D E\ .... 

f[A„=0]=( O D ) ^[An=l] = ( c pi U[An=-l]=( q ) • (4-8) 



The subscript involving An indicates the unique value of An which is involved in a history where 

'0 0^ 
v0 0, 



'0 0^ 

the initial and final states are (0-1,172) and (tTi,<7 2 ) respectively. = I _ _ I and A, B, C and D 



are 2x2 matrices. 
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4.1.2 Generating function method 

The evolution equation for V(rii,n2;t) is easily derived by considering the probability 
P(a- 1 ,(T2,n 1 ,n2;t) = ^ o 2 |U(ni, n 2 ; t)\a[, a 2 ) P(a[, o' 2 \ t = 0). 



(4.9) 



P(&l) 02) n\, ri2) t) is the probability that the system is in configuration (01,02) a time t and has 
received the heat amounts Q\ — — niAe and Q 2 — n 2 Ae during the time interval [0, t] when 
the initial probability distribution for the spins is P{a' 11 a' 2 ;t = 0). The evolution equation for 
P(o'i ) o'2,Tii,Ti2',t) is a generalization of the master equation (3.1l which governs the evolution of 
-P(ci, 02; t). By taking into account the explicit expression (2.9 1 for the transition rates we get 



2 dP{ax,a2,n-i,n 2 ;t) 

V\ + V2 dt 



- [1 - o-i<72(^i7i + ^272)] P(cri,a 2l n ll n 2 ;t) (4.10) 

-V x [1 + 0-1 0-271] P(-o 1 .cr 2 ,«i - oio 2 ,n 2 ;t) 
-V 2 [1 + 01O272] P{pu _cr 2, ni,n 2 + 01O2; t) 



where the dimensionless inverse time scales j/ a 's are defined in ( 3.4: 1 



The operator in the r.h.s. of the evolution equation (4.101 is partially diagonalized by con- 
sidering the generating function P(oi, a 2 , Zi, z 2 ; t) = J2n 1= -oo YZ^=-oo z i lz 2 2p { a h °2, "1, n 2 ; t) 
which is well defined for Z\ and z 2 with moduli equal to 1. Considering the latter generating 
function is equivalent to introducing 



+00 



+00 



d( Zl ,z 2 ;t)= ]T ]T *WLI(m,n2;t). 



(4.11) 



ni— — 00 n2——oo 



We notice that since P(ai,a 2 ,ni,n 2 ;t = 0) = S ni o§ n2 flP(ai, cr 2 ;t = 0), the initial value of 
V(zi,z 2 ;t) is merely \)(z\,z 2 ;t — 0) = D4, where D4 denotes the identity 4x4 matrix. The 
inversion formula which allows to retrieve V(ni,n 2 ;t) reads 



U(m,n 2 ;t) 



dzi 



ki|= 



dz 2 l)(zi,z 2 ;t) 



1 z7a J\z 2 \ = l z7rl z l z 2 



(4.12) 



The decomposition (4.6) of V(ni,n 2 ;t) leads to a similar decomposition for V(zi,z 2 ;t) 



V(zi,z 2 ;t)= Y ^lAn](zi,z 2 ;t). 
An=0,+1,-1 



(4.13) 



The decomposition ( |4.7| of U(ni, 77,2; i) into three 4x4 matrices (enforced by the constraint (4.5 1 
due to energy conservation) is also valid for V(zi,z 2 ;t). Moreover U[An] {z\ 1 z 2 \ t) has necessarily 
the following dependence upon z 2 and Z1Z2: U A"1 i z i, z 2 ; t) = z£ n V[A n ](ziZ 2 ;t) Therefore, by 



using the change of variable z\ — > z = z\z 2 in (4.121 one gets 
U(ni,m + An;t) - 



dz U [An] 0*1 = z,z 2 = l;i) 



>\ z \=i 2vri z™i+i 
4.1.3 Diagonalization 

The evolution of V(z\, z 2 ;t) with the dimensionless time variable 

v\ + v 2 



(4.14) 



reads 



dU(zi,z 2 ;T) 
dr 



= A(zi,z 2 )U(zi,z 2 ;t), 



(4.15) 
(4.16) 
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where, from the evolution equation (4.101 



{zi,z 2 ) = -h + 








b 2 







7* 


h 


h 


C2 


Cl 


~7* 







(-2 





-1*1 



(4.17) 



with the following notations: b\ = Ux(l + 7i)zi, b 2 = v 2 {l + 72)^1 c\ — ~~ 7i)^ anc ^ 
c 2 = V 2 (l - 72)22- 

Since the transition rates are invariant under the simultaneous changes of both spin signs, it 
is convenient to consider the transformed matrix 



with 



\'{ Zl ,z 2 ) = F- 1 A(z 1 ,z 2 )i 



/l 1 \ 

11 
1-10 

\o 1-1/ 



(4.18) 



(4.19) 



The change of basis corresponding to P 1 ensures that A'(z\,z 2 ) corresponds to two sets of de- 
coupled equations, 



K'(zi,z 2 ) = -I 



with the already defined matrix = ^ ^ and, for e = ±, 
trace of B, vanishes so that 



7* e6i + b 2 
eci + c 2 -7* 



]( Zl ,z 2 ) = A t ( Zl z 2 )i 



with 



where 



A e (z) = l-2A + e 



{A + B)z + (A — B)'- 



^ = ^1^2(1-7172) and B = viv 2 (71 — 72) . 
We notice that A > \B\. As a consequence, 

f e TB + (zi,z 2 ) Q 



e rA'(z u z 2 ) = e -r x 



3-(zi,z 2 ) 



where 



e rB e (z 1: z 2 ) = cosh ( TN /A e (^^))n 2 + 



h(r A /A 6 (z 1 g 2 )) [ 



<{zi,z 2 ). 



(4.20) 
. The 
(4.21) 

(4.22) 
(4.23) 

(4.24) 

(4.25) 



Moreover the eigenvalues of the matrix \{y\ + v 2 )k{z\, z 2 ) are, with the notations e = ± and 
V = ±, 



^\z) = V -^[-l + r,VMz) 



(4.26) 



4.1.4 Results for the generating function 



The relation (4.181 between A(zi,z 2 ) and A'(zi,z 2 ) allows to calculate V(z±, z 2 ;t) = e T ^ Zl > Z2 > as 



U(zi, z 2 \ t) = P 1 e rA ( z ^ z 2)p. From the explicit expressions (4.241 and (4.25) we get the matrices 



17 



U[An] z 2', t) defined by the analog of the matrix splitting (4.7 1. The 16 matrix elements can be 
written in the compact form 



Ol , (7 2 1 U [An=0] (zi ) 2 2 ; t) I a x , cr 2 ) 

(-01, -<7 2 |IU[A n =0] (*1 5 7") 1 0-2) 

(-CTi , cr 2 1 U [An ] {zi , z 2 ; r) I a x , cr 2 ) 



C + (ziz 2 ;t) + o-icr 2 7*5 + (ziz 2 ;r) 
C~(zi£ 2 ;t) + CTicr 2 7*S' _ (z 1 z 2 ;r) 
S„ 



(4.27) 



An=±l 



((7i, -cr 2 1 U [ A „] (zx,z 2 ;t) I ax, a 2 ) = S ai a 2 

An— ±1 



An zf n F+ n (zxz 2 ;r) 
An z£ n F7{zxZ 2 ;r) 



with 



C ± (z;r) = e- T 



1 r 



S ± (z;t) = e 



1 



sh(r^A + (z) ± cosh(r^A_(z) 
inh(Tv/A+(z)) sinh(rVA_(^)) 



± 



and 



^An( z i T ) = ^1 (1 " A "7i) 5 ± (z;r) + Z7 2 (1 - An 72 ) S^(z;r). 



(4.28) 
(4.29) 

(4.30) 



4.1.5 Results for the joint probability 



U(ni,ni + An;t) is derived from U[A n i {z\ , z 2 ; t) through the unit circle integral in (4.141. In fact, 
because of the parity property of the cosh and sinh functions, the functions C^(z; r) and S (z; r) 
are functions not of y/ A e (z) but only of A e (z). According to the definition (4.22) of A e (z), the 
single singular points of A e (z) are z = and z = +00, and the same is true for the integrands 
in § lzl=1 (dz/2iri)z-( m +VC ± (z;T) and f M=1 (dz/2iti )z~ ( - m+1 ^ S ± (z; r). The latter formulae involve 
integrals such as 

/ ttt — tt c °sh (r v/ATfTf) = f d6 

f M = 1 27ri 2™+! V V >) J_ v 



i m9 , 



2tt 



cosh tJ A_(e if? ) 



(4.31) 



By using the variable change 8 = 8' + ir and the periodicity of the integrand we get that 



dz 1 
As a consequence 

and 

where 

and 



cosh (rv/AT^)) =(-1)' 



=1 27ri z 



^COsh ( T yjA+{z) 



dz 

|*|=i 27ri z' 



j TI C ± (z;r)= 1±( o 1) '" e -c m (r) 



<iz 1 
|, ]=1 2^Tz^+i 

c m (r) = 



5 ± (z;r) = 



2 

i±(-iy 



e r s m (r), 



_i 27ri z' 



^cosh (rv/A+(z)) 



dz 1 sinh rv/A^ 



(4.32) 

(4.33) 
(4.34) 
(4.35) 

(4.36) 



Eventually, the matrix elements of V(rix,nx + An;t) are derived from the expressions ( |4.27| ) 
for the matrix elements of U[An] ( z 1j z 2! t) with the result 

(o-i, (7 2 1 U(n,n;t) I (7i, o- 2 ) = S cvcn (n) x C/ (0) (n, ctict 2 ; r) (4.37) 
(-cri,-cr 2 |U(n,n;i)|cri,cr 2 ) = J od d(n) x £/ (0) (n, (7icr 2 ; r) 
(-(7i,CT 2 |U(n, n + An; t) I £7i, a 2 ) = <5 CTl(T2 ,An x <J dd(n) x C/ (1) (n, An; r) 

An— ±1 

(cti, -CT 2 |U(n, n + An;i)|cri, ct 2 ) = J fflffJ ,A n x 5 ovcn (n) x (7 (1) (n, An; r), 

An— ±1 
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where 5 even (n) = \ [1 + (-1)"] and <5 odd (n) = \ [1 - (-1)"] while 

UW(n,a 1( T 2 ;T) = e~ r [c„(r) + ^a x o 2 s n {r)\ (4.38) 
[7 (1) (n, An;r) = e~ r [V 2 (1 - 72A71) s„(r) + V x (1 - 7! An) s n+An (r)] . 

We notice that the parity condition factors | [1 ± (—1)™] have a simple interpretation. During 
a history such that spin o\ is in the same state (in flipped states) in the initial and final con- 
figurations, the thermal bath 1 has flipped spin a an even (odd) number of times, so that the 
corresponding sum nAe of the successive amounts ±Ae dissipated towards the thermal bath 1 is 
necessarily an even (odd) multiple of Ae. 



4.2 Various explicit probabilities 



The probability that the system is in configuration (01, cr 2 ) at time t when the initial configurations 
are distributed according to the law Pq can be calculated, by virtue of the definition (4.111, as 

P Po (a 1 ,a 2 ;t)= ]T ^ (o-i,o 2 |U [An] (zi = M2 = l;t)\a[,a' 2 ) P Po (a[,a 2 ;t) . (4.39) 



The matrix elements of U[Ari]( z i = 1,^2 = are derived from (4.271 where, according to (4.221, 
A + (z = 1) = 1 while A_(z = 1) = 1 — 4A with A = V\V 2 (\ — 7172)- By using the identities 
Po(<ti,<j 2 ) + Pq(— a i, ~ a 2) = ij[l + criC27o], where 70 is the mean value of the product <ti<j 2 at 
time t = 0, 70 = (<Ji<J 2 )p , and P {ai,a 2 ) - Pq{-v\, -o 2 ) = \ [oi^o + 0- 2 (o 2 ) o ], where (cr'^o 
denotes the initial mean value of o\ , a straightforward calculation leads to 



Pp {o-i,(T 2 ;t) = ~[1 + 01027*] + ^cricr 2 [7o - 7*]e 2r 



(4.40) 



M^Do + 0-2(^)0] 



l + <ri0 . 2 ^) e -(l-«)r + ^ 



o"io-2 — e 

a 



(1+c*)t 



8 Q 



i/ 2 + CTlC2(^l7l - ^272)] 



e -(l-a)r _ p -(l+«)r 



where a = vl — 4 A. When Pq is invariant under the simultaneous reversal of the spins o\ and a 2 , 
Po(<Ti,a 2 ) = (1/4) [1 + <7i<j27o] and only the terms in the first line of (4.401 do contribute. Then 
the evolution of Pp cr 2 ; t) towards the stationary distribution P s t(c 1 , a 2 ) = ^ [1 -\- (?i& 2 7+ ] 
involves only one time scale, namely l/(v\ + v 2 ) (because r = \{y\ + v 2 )t). 

For any initial probability distribution Pq of the spins, the probability that at time t the system 
is in configuration (01,02) and has received a heat amount Q± = — niAe from the bath 1 and 
0.1 = n 2 Ae from the bath 2 is given by 



Ppo (0-1,0-2, rai, n 2 ;t) = ( a i^2\^(ni,n 2 ;t)\a[,(7 2 )Pp {a[,a 2 ;t) 



(4.41) 



From (4.371 we get that 

P Po (ai,a 2 ,n,n;t) 
Pp„ (01, o- 2 , n,n + An; t) 



An=±l 



J7 (0) (n,criCT2;-r) [S CVC n(n)P (ai, a 2 ) + S odli (n)P (-ai, -a 2 )} (4. 
<W 2 -AnP W (n, An; r) [<5 e ven(7i)-Po(0i, -cr 2 ) + <5 dd("0-Fo(~oi, 



Various joint probabilities can be derived from these expressions. 

The joint probability that at time t the system has received a heat amount Qi = — niAe 
from the bath 1 and a heat amount Q 2 = n 2 Ae from the bath 2 is calculated as Pp (711,712; t) = 
J2a u a 2 p Pa (o"i> a 2, m, n 2 ; t), with the result 



Pp (n, n; r) 
Pp (n, n + An; t) 



e T [c„(t) + 7* 70 s„(t)] 
1, 



(4.43) 



= e T -[l + An 70] [v 2 (1 - An 72) s n (r) + ^1 (1 - An 71) s„ +A „(r)] , 

An— ±1 z 
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where 70 has been defined before (4.401. 

The probability that at time t the system is in a configuration where <J\<Ji is equal to 1 and 
that the system has received a heat amount Q\ — — niAe from the bath 1 is defined as 

Pp (ata 2 =±l,ni;t) = ^2 S aia2> ±i ^ P Po (cti, cr 2 , «i, «i + An; i) (4.44) 

cri,(T2 An=— 1,0,+1 



From (4.421 we get that 



(4.45) 



Pp (cicr 2 = +l,ni; t) = 
V T {(1 + 70) c ni (r) + [2 7 * + (z/ 2 - vm) (1 - 70)] s ni (r) + Fi (1 + 71) (1 - 7o) V-lW) 



and 



Pp (cria2 = -l,ni;r)= (4.46) 
-e~ T {(1 - 7o ) c ni (r) + [-2 7 , + {V 2 + P171) (1 + 7o )] s ni (r) + V x (1 - 7i ) (1 + 7o ) s„ 1+ i(r)} . 

The expressions for Pp (<7icr2 = il,H2;r) are obtained from the latter equations by making the 
exchanges V\ O f72 and 71 -H- 72 and the replacements c ni — > c„ 2 , s ni — ¥ s„ 2 and s ni -i — > s„ 2+ i 
for <7i<72 — 1 an d s„ 1+ i — > s„ 2 _i for cricr2 = — 1- 

From these expressions we get the probability distribution for only one heat amount Q\ or Q2 



Pp {ni',r) 

and similarly 
Pp ( n 2-,T) 



e T {c ni (r) + [V 2 +V lllln ]s ni {T) (4.47) 
+ ^i(! - 7i)(l + 7o)s„ 1+ i(t) + -Fi(l + 7i)(l - 7o)s ni _i(r)| 

e~ T { c n 2 {T) + [Pi +P27270] s n2 (r) (4.48) 
+ 2^2(1 + 7a)(l - 7o)sn 2 +i(r) + *F 2 (1 - 7 2 )(1 + 7o)s„ 2 -i( r )}- 



Moreover, since J2n=-oo e Tc n( T ) = e T cosh A + (z = l)j = e r coshr, and similarly 

'Y^ti=-oo e ~ T s n{ T ) — e~ r sinhr, the probability that the total heat amount received from both 
thermostats is Q\ + Q2 — ("2 ~ «i)Ae reads 



Pp a ( n 2 - ni = 0; r) 
Pp ( n 2 -n 1 = An; r) 
As a consequence 



\ [1 + 7*7o] + ~ [1 - 7*7o] e 2r 
1 



(4.49) 



An=±l 4 



{1 - 7,70 + An [ 7o - 7,]} [1 - e- 2r ] 



([Qi + Q 2 ] T )p = g bo-7.] [1 



-2t 



] Ae, 



(4.50) 



and we retrieve property (3.44). 

We notice that all formulae are still valid in the limit where T\ vanishes, namely where /3iAe 
goes to infinity. 

4.3 Excess heats 

The excess heats associated with the transition between two different steady states have been 
defined in (3.481. For the two-spin system they can be explicitly calculated. Indeed, when the 
system is initially prepared in the stationary state with distribution Pq by contact with two 



20 



thermostats at the inverse temperatures (3® and /3° and when it is put at time t = in contact 
with two heat baths at inverse temperatures f3i and /3 2 , then the mean heat amount received from 
the thermostat 1 between time t = and time t is given (with convention (4.1 1) by 



(Qi(i)>P =-Ae "i E E E ( CT i' CT 2|U[A„](m;i)ki,a 2 ) 



(4.51) 



rii=— oo An=0,l, — 1 fi ,iT2 cr! ,cr» 



By virtue of the definition (4.11 ) of the relevant generating function and the decomposition (4.131, 
the latter expression can be rewritten as 

d 



(fix(*))j 



-Ae 



E EE { a 'i^'-iP[An}{zi,zi]t)\a 1 ,a2) 

An=0,l,— 1 cri,CT 2 <r;,er 2 



{Z! =l,Z2 = 1). 



(4.52) 

From the explicit expressions (4.271 for the matrix elements and the value (4.221 of A + (z), a 
straightforward calculation leads to 



(Qi(t)) Po = (ji)rt* + [TO (7i - 7 2 ) ~ Fi (71 - 7*)] 



Ae 



X _ e -(^l+t>2)t 



(4.53) 



where (ji) a t — — 0*2)st is given in (3.271 and 7° is a function of i/i, v%, ft® and written in (3.51. 
A similar calculation yields 



(Q2(t))p = 0*2>st* + [-V1V2 (71 - 72) - V 2 (72 - 7*)] 



Ae 



I _ e -(fi+fa)t 



With the sign convention of definition (3.481 



T lift, ft] 

-) I [ft .ft] 
« exc - 2 l[^,/3g] 



[j/ x i/ 2 (71 - 72) - v 1 ( 7l - 7°)] 



Ae 



(4.54) 



(4.55) 



[-V1V2 (71 - 72) ~ ^2 (72 - 7*)] 



Ae 



Comparison with the expressions (3.141 for the stationary mean energies in the initial and final 
stationary states shows that the excess heats coming from both thermostats indeed satisfy the 
identity (13.49 1. 



4.4 Symmetry property for reversed heat transfers (when T\ 7^ 0) spe- 
cific to the model 

4.4.1 Symmetry arising from the modified detailed balance 

The symmetry properties for reversed heat transfers when T\ ^ are more conveniently exhibited 
after a change of variable in the complex plane where the integrals involved in Pp (Qi, Q.2,t) are 
defined. The functions at stake c„(r) and s n (r) are defined in (4.351 and (4.36) while A + (z) is 
given in (4.221. The origin z = is a singular point in A+(z) and zA + (z) = (A + B)P(z) where 
the second order polynomial P(z) — (z — z + )(z — zJ) vanishes for two roots z + and z_. The 
product of the roots is equal to 



A — B 
A + B 



with 



p : 



3 (/3 1 -/3 2 )Ae/2 



(4.56) 



and the sum of the roots is equal to — (1 — 2A)/(A + B). 

When T\ ^ 0, z + z_ does not vanish and by using the variable change £ 



z I -JZ-Z+, namely 



C = zp, 



(4.57) 
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the unit circle is changed into a circle with radius p, while the roots z + and z_ are changed into £_ 
and £ + with C-C+ = 1- Then, for a function such as cosh(-y/A + (z)) or sinh(r-y/ A + (z))/ yj A + (z), 
each of which is in fact a function of A+(z) denoted by f(A + (z)), 



' K ^r/(A + (C)) 



with A_|_(C) = A + (£/p). A + (£) is a symmetric function of £ and l/£ 

A + (()=b + a 1 ^- 

with a = 2\ / A 2 - B 2 and 6 = 1 - 2A, namely 

a = 2F 1 F 2V /(l- 7l 2 )(l- 72 2 ) 
6=1- 2i7 1 i7 2 (1 - 7172) • 



(4.58) 



(4.59) 



(4.60) 
(4.61) 



circle |C| = P can be deformed into the unit circle and we get the identity 



Since the only singular points in the integrand in the r.h.s. of (|4.58| are ( = and £ = 00, the 

(4.62) 



| z |=i 2tti z" +i J| Z | =1 2tti z n + L 



By inserting the latter identity in the definitions (4.351 and (4.361 we get the relations 



c„(r) 
s n (r) 



P n c n {r) 
P n s n {r) 



where 



and 



c„(t) = 



s n (r) 



dz 1 



=1 27ri z n+l 



--f cosh ryA+(z) 



1 sinhlryA + (z 



Since A + (z) is invariant under the exchange of £ and l/£, 

c„(r) = cj n |(r) 
s n (r) = S|„|(r). 



(4.63) 



(4.64) 



(4.65) 



(4.66) 



Therefore the functions involved in the matrix elements (4.371 of V(nx,ri2\t) can be rewritten as 

U (0) '(n,o-i<7 a ;7-) = e~ T p" [c w (r) + 7*cxi<7 2 sj„| (t)] (4.67) 
C/W(n,An;r) = e" T p" [z7 2 (1 - 72 An) s w (t) + p t (1 - 7l An) p a "s|„+a«|(t)] . 

where p. defined in (4.56), also reads p = \J{\ +7i)(l — 72)/(l — 7i)(l +72)- 

According to paper I , the modified detailed balance implies some time-reversal symmetry 
property for histories, which itself entails some relation between probabilities of forward and 
backward evolutions where the given initial and final configurations are exchanged (and the heat 
amounts are changed into their opposite values). In the spin model language, with the definitions 
(4.1), the symmetry exhibited in paper I for the probability that the system evolves from an initial 
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configuration Co = (ai, <7 2 ) to a final configuration = (o^, "2) while receiving the heat amount 
Qi = — niAe and Q 2 — «2^ e reads, for non-vanishing matrix elements, 



(a' 1 ,a 2 \l)(n 1 ,n2;t)\a 1 ,a 2 ) 
(ci, cr 2 |lU(-ni, -n 2 ;t)\cr' 1 ,a 2 ) 



Comparison of the latter relation with the expressions (4.371 implies that 

U-C0)( n ,tri(72;r) 



t/( )(-n 7( T 1( 7 2 ;T) 



,[n/9i-n^ 3 ]Ae 



and 



= e [nft-(n+A f 



(4.68) 



(4.69) 



(4.70) 



The latter relations can be checked from the explicit expressions (4.67). 

We notice that the relation ( |4.68[ ) can also be retrieved by noticing that the modified detailed 
balance entails that the matrix A(z\,z 2 ), which rules the evolution of V(zi, z 2 ;t) according to 
(4.161, obeys the symmetry k(z\, z 2 ) — A T (e~^ lAe / z\, e^ Ae / z 2 ), where A T denotes the transposed 
matrix of A. Therefore, after the variable change Z\ — Ci/pi an d z 2 — C2/P2 with pi — exp(/3i Ae/2) 
and p 2 — exp(— /3 2 Ae/2), the matrix A(zi,z 2 ) becomes the matrix A(C 1 ,C 2 ) = ^{ z i — Ci//°i; z 2 — 
C2/P2) which obeys the symmetry A T (d, £ 2 ) — A(l/£i, l/Cz)- Then the derivation of the symmetry 
(4.681 is the following. First we make the variable change z\ — (i/pi and z 2 — C2/P2 in the integral 
representation (4.121 for V(ni,n 2 ;t). Since A(d,£ 2 ) has no non-analyticity, apart from the l/d 



and I/C2 singular terms, the integrals on the circles with radii equal to p\ and p 2 are equal to the 
integrals with the same integrands on the circles with radii equal to 1. If C is on the unit circle, 
1 /£ is also on this circle and we can make the variable change d — > and C2 — > I/C2 ; then the 
symmetry of A leads to the symmetry (4.68). 



4.4.2 Symmetry specific to the model 

The present model happens to obey a very specific symmetry property for reversed heat transfers 
when the initial state of the system has an arbitrary probability distribution Pq. According to 



(4.371, after summation over the final configuration, 
^2 Wx,^2p{ n ,n;t)\ax,(72) 



{<J / 1 ,a / 2 \V{n 7 n + An;t)\a 1 ,a 2 ) 



An=±l 



(T2 , An ■ 



(4.71) 



Therefore, when the initial configurations are distributed with an arbitrary probability Pq 

^K,a 2 |U(n,n;i)|P ) = Yl ^ K^icr 2 ; r)P (tn, ct 2 ) (4.72) 



and 



V (a' 1 ,a' 2 \V(n,n + An;t)\P a ) = U {1) (n, An; r) x P 

(ci<72 = An) (4.73) 

z — ' An— ±1 



where P (<7i<t 2 — An) denotes the probability that the product CTi<7 2 is equal to An in the initial 
configuration. 



Then from the expressions (4.721 and (4.731 for Pp a (Qi, Q,2,t) in the cases Qi — —Q2 and 
Qi ^ —Q.2 respectively, and by virtue of the consequences ( 4.69| -(4.70 1 of the modified detailed 
balance, we get the property 



Pp (Qi,Q2,;t) 
Pp {-Qi,-Q2-,t) 



_ e -/3iCi-/3 2 C2 



^Qi+Q 2 ,0 + X! 



+ Q 2 ,cAe 



€ = ±1 



P{oiQ2 = e;t = 0) 
l-P(cr 1 (T 2 = e;i = 0) 



(4.74) 
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We stress that this relation is very specific to the present model. In the cases where the initial 
distribution is equal either to an equilibrium state distribution at the inverse temperature /3q or 
to the stationary state distribution which is a canonical distribution at the inverse temperature 
/3*, the relation allows to retrieve the generic relation studied in paper I. 

4.5 Decaying property of joint probabilities for large heat exchanges 

All quantities of interest involve the coefficients c m (r) and/or s m (r), computed via contour in- 
tegrals in ( |4.35[ ) (4.361. The integrands involve functions which are holomorphic for z g C*, the 



pointed complex plane: the quantity A + (z) defined in (4.221 has this property and the func- 
tions cosh(u;) and smh(w)/w are entire even functions, so that the square roots in the composed 

functions cosh A + (z)j and sinh (t^ A + (z)^j / yj A + (z) do no harm. 



Hence, in the formula? for c m (r) and s m (r), contours can be deformed. For r e]0, +oo[, let 
C(t, r) = sup| z | =r | cosh (t^/ A + (z)j \ < +oo and S(t, r) = sup| z | =r sinh (t^J A + (z)j / y/ A + {z) 
+oo. 

Taking \z\ = r as integration contour, one gets immediately that, for each r, c m (r) < 
C(T,r)r~ m and s m (r) < S(r,r)r~ m . This shows that c m (r) and s m (r) are o(e _if ' m ') at large |m| 
for any K . 

With some efforts, we could get some explicit upper bounds for C(r, r) and S{r,r). Then we 
could extremize over r to get a subexponential bound for c m (r) and s m (r), but we shall not need 
this refinement. 

In the limit Ti = 0, A + (z) is in fact holomorphic for z € C so that c m (T) and s m (r) vanish for 
m= -1,-2,- 



< 



5 Heat amount cumulants for any T\ and 

5.1 Generic properties for a system with a finite number of configura- 
tions 

5.1.1 Characteristic function for the heat amount Q2 

The random variable Q 2 can take only discrete values n 2 Ae, where n 2 is a positive or negative 
integer. Therefore its probability density n(Q; t), defined as II(Q; t)dQ = P (Q 2 G [Q, Q + dQ[; t), 
reads 

1 +°° / q \ 
n ( Q ;*)=A^ E s(^--n 2 \P(n 2 ;t), (5.1) 

ri2 — — 00 ^ ' 

where S stands for the Dirac distribution. Since P(n 2 ; t) decays faster than exp(— i4T|n 2 |) for any 
K > when |n 2 | goes to infinity (see subsection 4.5 1, the Laplace transform G(A;i) of II(Q), i.e. 
the characteristic function of the random variable Q 2 , is well defined for any A, 

/+00 
dQe XQ U(Q;t) = (e XQ ^). (5.2) 
-co 

As a consequence, the properties of the probability density H(Q;t) can be investigated through 
those of its Laplace transform, thanks to the inversion formula 

n(Q;i)=/ £± e - xa G(\;t) (5.3) 



J X 



2?ri 



According to the property (5.1), G(A) is a periodic function of A with period equal to i27r/Ae, 
and the r. h. s. of the latter formula can be written as 



1 d9 



e -iea/AeQ(^_ A y e -i2, mQ /Ae (5 4) 

Ae./n 2tt Ue / ^ 



m— — 00 
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By virtue of Poisson equality X)m= 
leads to 



P(n 2 ;t) 



i27rmQ/Ae = J^t^L-oc & - "2), comparison with (JKTl 
2 " dO 



2tt 



2 G 



Ae' 



(5.5) 



The latter equality coincides with the inverse formula (analogous to (4.12)) in terms of the gener- 



ating function, G(z;t) — J2n°^ 



2 P(n2-,t), 



P{n 2 -t) 



dz 



1 



=1 27ri z n * 



(5.6) 



where G{z;t) = G((l/Ae) Inz; t). 



5.1.2 Relation between long-time cumulants per unit time for Q\ and Q2 

The generic properties of the cumulants of Q± and Q2 have been reviewed in paper I . We recall 
those which will be useful in the following. For a Markov process, the long-time behaviors of these 
cumulants are proportional to the time t elapsed from the beginning of the measurements. The 
asymptotic behavior of the cumulants per unit time are given by the derivatives of 



«a(A) 



lim -ln(e XQ "^) 

t— >+oo t 



according to 



>] 

lim 

t-¥ + oo t 



d p a a {\ a ) 



dX P a 



for a = {1,2}. 



(5.7) 



(5.8) 



A a =0 



Moreover, in the case of a system with a finite number of configurations, Q\ + Q2 = £(Cf) — 
f (Cq) is restricted to some finite interval and 



ai(A) = a 2 (-X). 

As a consequence the long-time cumulants per unit time obey the following relations 



M 



lim 



M 
"Q 2 



(-lriim 

t— too t 



t^roo t 

5.2 Explicit formulae for the cumulants per unit time 



(5.9) 



(5.10) 



According to the relation (5.10 1 between the long-time cumulants per unit time for Qi Q 2 , we have 
only to consider the cumulants for the heat amount Q2 received from the bath 2. For the two-spin 
system, where Q2 = «2Ae, it is convenient to introduce the cumulants k„1 for the dimensionless 
variable ri2 and the associated characteristic function (e A ™ 2 ), where A is a dimensionless variable. 
According to the relation (5.8) the long-time behavior of the cumulants per unit time are derived 
through the relation 



1 



lim - k_ 



Ip] - 



dPa 2 (X) 



x=o 



with 



012(A) = lim - 
t^+00 t 



ln(e A " 2 ) 



(5.11) 



(5.12) 



According to the definition of (e A ™ 2 ), the relation (4.9 1 between the probability P{<J\, 02, ni, n?; t) 



and the operator U(ni,n2;i), together with the definition (4.111 of V(zi,Z2',t) and its evolution 
iracterist 

Jn2 ) = E E 



equation (4.161, the characteristic function may be expressed as 

I 0-1,0-2 1 e 



5-tA(zi = l,z 2 =e A )| 



(5.13) 
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According to (4.171 



(zi = l,z 2 = 6 



(1 + 7*)D 4 is a real positive matrix and the Perron- 



Frobenius theorem can be applied. Henceforth 02(A) coincides with the eigenvalue of the matrix 
\(h>i + j/2)/2]A(zi = 1,2:2 = e A ) with the largest modulus (and which is necessarily real). The 



four eigenvalues to consider are the n^'^(z = e A )'s which are given by the expression ( |4.26[ ), with 
e = ± and r) = ±. The one with the largest modulus corresponds to (e, if) = (+, +) and reads 



a 2 (A) = 



-1 



1 - 2A + (A + B)e x + (A — B)e~ 



(5.14) 



where A and B are defined in (4.231 



It is plain to calculate a number of cumulants from (5.111. Their behavior as a function of the 



kinetic parameter ^2/^1 exhibits some interesting features. For large ^2/^1, the cumulants go to 
a limit which is the same for all odd and for all even cumulants. as will be explained in section 
[7] Fig[l] illustrates this convergence, which gets slower and slower for higher moments. The 
first six cumulants are represented. This figure also shows some oscillations at finite ^2/^1- These 
oscillations become more and more visible on higher cumulants. Fig|2]illustrates this phenomenon. 
Cumulants from the fifth to the ninth are represented. In both figures, the other model parameters 
are fixed to the sample values 71 = 0.7, 72 = 0.4, v\ = 2. 
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Figure 1: First cumulants as a function of ^2/^1- Illustration of the large z^/^i behavior. The 
other model parameters are fixed to the sample values 71 = 0.7, 72 = 0.4, v\ = 2 



Only the first cumulants have analytic expressions simple enough to fit on a line. For the sake 



2G 




Figure 2: Higher cumulants. Illustration of the oscillations at finite Vi/v-y. The other model 
parameters are fixed to the sample values 71 = 0.7, 72 = 0.4, v\ = 1 



of conciseness, the results are first expressed in terms of the dimensionless time r = \(y\+ v<i)t as 

1 



lim -kW 

t->-+oo r 

lim i/eP] 
lim -iS 
lim — 



I? 

A-B 2 

B [l — 3^4 + 3B 2 ] 
A - 3A 2 + B 2 [-4 + 18A - 15B 2 ] 



(5.15) 



r— f+oo 7" 

All odd cumulants are proportional to B, because all odd powers of A in the expansion of the 
expression 



5.14) for (22(A) are proportional to B. The first three cumulants are rewritten in 



terms of the model parameters as 

(Q 2 )st 



lim 



lim 

t— >+oo t 

(Qi)st - (Q 2 ) 2 t 



hm 

t— >+oo £ 



- - / A"i + ^)Ae . 

f 1^2 (71 - 72) ^ = 2 )st 



2 , ( I / 1 +^)(Ae) 2 



^1^2 [l - 7i72 ~ ^1^2(71 - 72) 



^1^2 (71 - 72) 1 - 3i^i^2(l - 7172) + 3^2 (71 - 72) 



(5.16) 
(1/1 + i/ 2 )(Ae) 3 



(Q2) c i s the third cumulant, which is equal to the third centered moment, namely (Q| 

([Q2-(Q 2 )] 3 ). 



3\c — 
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At equilibrium 71 =72, so that B = 0: then, by virtue of the remark after (5.151, the long-time 



behavior of all odd cumulants of Q 2 is subdominant with respect to the elapsed time t. and in the 
long-time limit P(Q.2',t) becomes an even function of Q 2 at leading order in time t. The fourth 
cumulant of the cumulated heat Q2 received from the thermostat 2 per unit time does not vanish: 
limt_j._i_.oo j ln(e A ® 2 ) eq is not quadratic in A, and even in the long time limit the variable Q2 has a 
non-Gaussian distribution, contrarily to the variable [Q 2 — (22)] /y/i (for which all cumulants of 
order larger than 3 vanish in the infinite time limit). The first two even cumulants per unit time 
read 



lim 

t— . +00 



\£_)cq 



Vl v 2 (l - 7 2 ) 



(fi+" 2 )(Ae) i 



(5.17) 



lim 



/Q 4\c 



2/ cq 



= ^ 2 (l-7 2 ) [1- 37^ 2 (1 - 7 2 )] 



(v 1 +is 2 )(Ae) 4 



\Q\) C denotes the fourth cumulant, which can be expressed as (Qj) c = ([Q 2 - (Q 2 )] 4 )-3([Q 2 
For a system weakly out of equilibrium the Green-Kubo relation, namely 



(/?: 



,/. 3 )-. (/.,/.) 0i - ft 



- lim 

2 t— >-+oo 



\2 



2 \ 

2/ cq 



«Qa 



7 1 



(5.18) 



is indeed obeyed by the system, as it should be. This can be checked by comparing the expression 
(5.171 with the limit obtained when (/3i,/3 2 ) — > (0,0) for the ratio (j 2 )st/(0i — 02) which, by virtue 



of(3.27l, reads 



(5.19) 



(j2)st 71-72 (^1 + ^)(Ae) 2 

— V^Uo . 

Pi-fo Pi -fa 2 

When the system is far from equilibrium, comparison of the latter expression for (j2)st/(/3i — 02) 
with the expression (5.161 for the long-time limit of the second cumulant per unit time shows that 
0_)st/(/3i ~02) ^ limt-H-oo [( Qpgt - (Qj^st.A 1 ^ ^ should be (see subsection 5.3 of_ paper I). 
Indeed, by virtue of equation (5.141, a 2 (A) obeys the symmetry relation a 2 (A) = a 2 (— J 7 — A) with 
T = \n(A + B)/(A — B) = {j3\ — /3 2 )Ae but a 2 (A) is not a quadratic function of A, i.e Q 2 has a 
non-Gaussian distribution in the long-time limit. 



6 Large deviation function for the cumulated heat current 

Q2A 

In this section, we derive the large deviation function for the cumulated heat current Q 2 /£ by three 
methods. The first one is based on the general theory of large deviations for the definition of large 
deviation functions and uses one of its cornerstones, the Gartner-Ellis theorem. The second and 
the third rely on the fact that Q 2 takes discrete values in a t-independent set, and uses an ad-hoc 
definition of large deviation functions (see Appendix D of paper I ) . Though the general theory of 
large deviations and the ad-hoc definition for discrete exchanged quantities do not have to be the 
same, the ad-hoc definition is nevertheless a sensible definition of large deviations. Physically, the 
general and the ad-hoc definition are expected to yield the same result in a case as simple as the 
two-spin system, and our explicit computations can be seen as a proof of this fact. A natural tool 
to compute the ad-hoc large deviation function is via a contour integral representation, but as we 
shall see below, this method is surprisingly tricky even for the simple two-spin system at hand. In 
contrast with the general theory of large deviations, the contour integral method is the basis of 
a systematic expansion at large times. However, corrections are less universal than the dominant 
term. 

The cumulated heat current received from the heat bath 2 during the time interval t takes 
the values J = <2 2 /i, with Q 2 = nAe, n integer. By dimensional analysis, the large deviation 
function /g, (_7), which has the dimension of an inverse time, must be a function of 
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and we shall often consider the expressions of 



fa 2 (j) = fa 2 (J) (6-2) 

rather than those of /g 2 (i7) — f q 2 (J / /S.e) . Moreover, the explicit calculations are more conve- 
niently dealt with if, instead of considering j, we introduce the dimensionless current j associated 
with the dimensionless time r = [(z^i + ^ 2 )/2]i, 

J=- = r^—J- (6-3) 

The dimensionless large deviation function /g of j is such that tfQ 2 (j) = t/q (j), and the 
expression of /q 2 (j) can be retrieved from that for f q 2 (j) through 

/ e ,0> = ^/ ai (^j. (6.4) 

We notice that large deviation functions for other cumulated quantities are related to /g 2 . 
Indeed, in a system with a finite number of configurations Qi + Q 2 is bounded and, as a consequence 
of the general theory of large deviations (see e.g. paper I), 

faAJ) = fa 2 (-J)- (6-5) 

In the same vein, as A cxch S = P1Q1 + P2Q2 = —{Pi ~ fii)Q,2 + fii(Qi + Q2), the large deviation 
function for A^chS" and that for Q2 satisfy the simple relation 

fA^s(J) = f Q2 {-^ 2 ). (6.6) 

6.1 Derivation from Gartner-Ellis theorem 
6.1.1 Method 



By analogy with (5.121, we introduce the dimensionless function 

1 • ' >2 



a 2 (A) = lim -ln(e A " 2 >. (6.7) 



A simplified version of the Gartner-Ellis theorem (see e.g. the review for physicists [5T] or the 
mathematical point of view [22]) states that, if «2(A) exists and is differentiable for all A in IR, 
then the large deviation function of the current j = 112 /t exists and it can be calculated as the 
Legendre-Fenchel transform of «2(A), namely, with the signs chosen in the definitions used in the 
present paper, 

7 Q2 (j) = mm{a 2 (A)-Aj}. (6.8) 

AgR 

As a consequence, since a 2 (A) obeys the symmetry a 2 (A) = 5 2 (— (/3i — /3 2 )Ae — A) (as can be 



checked from ( 5. 14 1) , f(j) obeys the fluctuation relation f(j) — f(—j) = — /3 2 )Aex j. Moreover, 
the cumulant generating function (e A ™ 2 ) is necessarily convex (downward). More precisely, in the 
present case a 2 (A) is not only continuously differentiable but also stricly convex for all real A, so 
that the minimum in the definition of the Legendre-Fenchel transform can be readily calculated 
by using the Legendre transform, 

7q 2 (J)=32(A c (j))-jA c C?) with %A C )=J. (6.9) 
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6.1.2 Various expressions for Jq 2 and its properties 



From the relation 02(A) = [{y\ + i^/^Ja^A) and the expression (5.141 for a 2 (\), when T\ 7^ 
(71 7^ 1), A 7^ B and we get 



VA 2 - B 2 



-1 + v/i-2A + 2r(j). 



(6.10) 



cosh 1 a; denotes the positive real whose hyperbolic cosine is equal to x. namely cosh 1 x — 
In [x + \Jx 2 — \\ , and 

y(j) - j 2 + ^ + (1-2^ + ^2-52. (6.11) 



The expression for /g (j) involves the combinations of the model parameters 

In^l^f =(^1-^2)^, A 2 -B 2 =I72l72(l- 7 2)(l_ 7 2) j A = P 1 P 2 (l- 7 nfc). (6.12) 



The expression (6.10 1 for /q 2 (j) can be rewritten in two different forms according to the sign 
of j. By using the identity cosh -1 x = In [x + \/x 2 — lj , the cosh -1 term in (6.10 1 can be split 
into two contributions and, according to the sign of j, we get 

IqM = -JHA - B)+jbx W(j) + ^Y 2 (j)-(A 2 -B 2 )} - 1 + y/l-2A + 2Y(j), (6.13) 



3<0 

while 

TqM = +jHA + B)-j]n 

3>0 



Y(j) + ^Y 2 (j)-(A 2 -B 2 )} - 1 + ^l-2A + 2Y(j). (6.14) 



In the limit where T\ vanishes (A — > B), the latter expressions yield the results discussed in section 

m 

The thermodynamical and kinetic parameters of the heat baths are disentangled if, in place of 
A and B, we consider the parameters 



p + = -(1 + 7i)(l - 72) and p_ = -(1 - 7i)(l + 72). 



(6.15) 



2" ' ' v ' ' 1 2 

The relations with j4 and -B are A = V\v 2 (p+ + p_) and £> = (p+ — P-)- Therefore, 



+ = p+/p_, V^4 2 - = 2i7 1 F 2 . v /p+pr. Then, by virtue of the relation ( |6~4 l 

-1 + v/1 - 2FiT7 2 [p+ +P-- 2Z(j)] . 



and the expression (6.10) for /q (j), /a 2 (j) reads 



/ C2 (j)=iW— -lilcosh- 1 



V\ + V 2 



where, with the definition Y{j)j\J A 2 — B 2 = Z(j)/ ^/p + p_, 
1 



2/ + v/4j 4 + [(^1 + v 2 ) 2 - 2v 1 v 2 {p++p-)]j 2 + (z^i^) 1 



P+P- 



(6.16) 



(6.17) 



By virtue of the definitions (6.15) of p + and p_, the thermodynamic parameters of the thermal 
baths appear in /g 2 (j) through the following combinations 

y|± = e ^-«^, VP^=^,\/(l-7i 2 )(l- 72 2 ), p++p- =1-7172- (6.18) 

At equlibrium the large deviation function is even. As 71 — 72 increases, the large deviation 
function becomes more and more asymmetric. In the zero temperature limit 71 = 1, the large 
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deviation function becomes infinite for j < 0. Fig[3] illustrates the changes in the shape of the 
large deviation function, with increasing departure from equilibrium. 

Some generic properties of a large deviation function can be checked in the case of the above 



explicit formulae. By virtue of ( 5.15| ) (j) st = B and one checks that /g 2 ((j) s t) = 0, /g 2 ((j) s t) = 
and 7g 2 (0) st ) = -1/(A - B 2 ), namely 



/o 2 ((j)st) 



lim 

r— >+oo T 



(6.19) 



fi-i)^f- and an even function of j. As 



The expression (6.101 for /g 2 (j) is the sum of a term j{(3\ 
a consequence, we check again that /q 2 (J) — /q 2 ( — j) = — /?2)Ae namely, by virtue of (6.4 1, 
Jq 2 {J) obeys the fluctuation relation /q 2 (i7) — /q 2 (— J) — {Pi — fti)3 ■ Moreover the absolute 
value of j in the expression (6.10) for /g 2 (j) is responsible for a (rather mild) singularity in the 
curve /q 2 (i7) at J = 0: a jump in the third derivative. 

We notice that the large current behavior of /g 2 (j) in the present model reads 



IJI- 



"2|j| In |j|. 



(6.20) 




Figure 3: Effect of increasing the non-equilibrium driving parameter 71—72- A few large deviation 
functions are represented. The leftmost curve is the equilibrium large deviation function and the 
dashed curve is the quadratic with the same curvature at the origin. The rightmost curve is for 

71 = 1 i.e. heat bath 1 at zero temperature. The other parameters are fixed to the sample values 

72 = 0.1, Vi = 100, Vx = 100. 
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6.2 Derivation from a saddle-point method 

Before embarking on the derivation, let us explain why the saddle point method is not straight- 
forward for this model. 

The saddle point approximation or expansion is well-suited for the asymptotic study of integrals 
of the form J r dz ty(z)e T ®W where r is some large real parameter and the functions ^(z), &(z) are 
holomorphic in a domain large enough that the initial integration contour T can be deformed to a 
steepest descent path while remaining within the holomorphicity domain during the deformation. 
One may also have to encircle some singularities when deforming the contour, and then one must 
keep track of their contributions, which may or may not dominate the saddle point contribution. 
This can of course be generalized to a finite sum J r dz 1 ^ i (z)e T ' Pi ^ when each individual term 
satisfies the hypotheses above. Note however that to get the leading behavior one may have to 
take into account possible destructive interferences between different pieces, for instance if the real 
parts of saddle point values are the same for several ^Vs, or if the saddle points for certain terms 
compete with encircled singularities for other terms. 

In our case, we deal with an integral of the type J r dz i ^i(z)e T ^ i ^ z ' where there is a single 
integration contour, and the su 

nf]Ei*i(^)e r * iW has nice holomorphicity properties that allow 
to deform contours (almost) freely, but each term in itself has singularities and cuts. So we have 
to face a kind of dilemma: either we want to keep holomorphicity, then the large parameter 
does not appear in an exponential - and to our knowledge no straightforward constant phase 
technique applies - or we look at each pure exponential piece individually, and then some branch 
cuts may prevent from deforming the contour purely as a constant phase steepest descent path: 
the steepest descent path is not closed, some parts of the original path are deformed along the 
cuts and they may dominate the saddle point. But also, the contribution of the pure exponential 
pieces may interfere. In our case, we have managed to show that in fact the interferences between 
contributions of one pure exponential and cut contributions from another pure exponential are 
destructive (with reminder terms controlled explicitly), leaving the contribution of only a single 
saddle point (not one saddle point for each pure exponential). But our argument relies on some 
tricks and features that appear to us at this stage as coincidences: we have not been able to identify 
a general framework avoiding our tedious analysis. And indeed, examples are known |23ll2"l] where 
(depending possibly on parameters) the cut contributions do or do not dominate the saddle point. 

To conclude these comments, let us mention one general direction that seems worth pursuing, 
though we have not been able to use it to simplify significantly our argument even in our simple 
case. In physical problems, the functions $i(z) will often be closely related to the different branches 
of a single algebraic function. For instance, the functions $i(z) are often closely related to the 
eigenvalues of some z-dependent matrix. So a natural route would be to regard the integrals not 
in the z plane, but on the approriate uniformizing Riemann surface, in our case an elliptic curve. 

We now turn to the detailed analysis. 

6.2.1 Method 

The current probability density ±I(J;t) is related to the probability P(Q2/(rAe) e + <£?[;t) 
by the definition P(Q 2 /(TAe) 6 \j,J + dj[;T) = H(j;T)dj. Since Q 2 /Ae can take only integer 
values, its density distribution II(j;t) is a sum of Dirac distributions 

n(j;r)= g S (j-^)rp(^ e =n;r). (6.21) 

In the long-time limit 

n(j;r) ~ Yl ^(j~-)^0^) (6-22) 



2 Which in our case consists of only two terms. 
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where P as (j;r) is a function of the continuous parameter j that we shall compute below, and 
which is such that the following asymptotic behavior holds: 



_ . r— >-+oo 
rj integer 



P as (j;r). (6.23) 



The notation g(j, t)\ t - intGger is a reminder of the rule that if the function g(j,r) is given by an 
integral representation, the latter must be calculated in the case where tj is an integer. By using 
one of the ad-hoc definitions of the large deviation function introduced in Appendix D.2 of paper 
I, the function P as (j;r) can be rewritten as 

P as (j;r) = A(j,T)e r7 ^ with lim - In r) = 0. (6.24) 

When one is interested only in the large deviation function, the only information to be retained 
from the latter equation is merely 

7q 2 © = T ^ O0 \ ln P (Q2 = rjAe; r)| rJ integcl . . (6.25) 

Consequently, /g (j) can be investigated by means of a saddle-point method applied to the 
representation of P{Q,2 — TjAe; t)\ t - intcgor in the complex z plane given by ( |4.48| ). In the latter 
expression P{Q,2 = nAe; r) is equal to e~ r times a linear combination of c„(r), s„(r), s„+i(r) 
and s„_i(r). When T\ = the expressions ( |4.35 1 and (4.361 of the latter functions are convenient 
for studying the large t behavior of c n=T j(r), s n=T j(r). When Ti 7^ the study is slightly more 
complicated and it is more conveniently performed by considering the related coefficients defined 
by c„(r) = p n c n (T) and s„(t) = p n s n {T) where the expression (4.561 of p is finite when T\ =^ 0. 
We present the details in the case where T\ ^ 0. 

When T\ 7^ 0, in the long-time limit, we have to consider the behaviors of the functions 



K c (j;t)=^ t j(t) = f ^A-l-coBhl r V 'A_C) ) ((>.2(i) 

and 

„ , sinh I t^/A+(z) 

^ 6 (An) (J; r) = ^ + a„ (r) = i £ ^ V (6.27) 

with An = 0, 1, — 1. It is sufficient to exhibit the derivation of the long-time behavior of K c (j; r), 
because the calculation of the long-time behavior of i4fi A ™- > (j; r) follows the same lines. Moreover, 
according to the property c„(r) = ci n i(r), we have to consider only the case where j > 0. 

For the study of the large r limit, the cosh function in the integrand of K c (j; r) is split into 
two exponentials, and K c (j;t) appears as the sum of two integrals 

2K c (j; r) = / J^ e r^) + i ^_ e r.W (.) (6 . 28) 

where 




$( ± )(z) = -Jlnz±W6+- z+-l. (6.29) 



We notice that, since tj is in fact an integer, exp(— rjlnz) is single valued and there is no cut in 
the complex plane z associated with the logarithmic function. However, since the cosh function 

has been split into two exponentials, we have to consider the two cuts associated with y A + (z). 
These cuts are 

] — 00, — x>] and [— ie<,0] (6.30) 

where — x > and — x < are the two negative real roots of the second-order polynomial zA + (z) where 
A + (z) is given in (4.591. The roots are such that < £< < 1 < x>. 



33 



6.2.2 Deformation of contours 



The large r behavior of K c (j; r) can be investigated by applying the saddle-point method to the 
contribution from the integral involving ^ + \z). For that purpose we have to find a way to deform 
the unit circle into a contour that goes through a saddle point along a constant phase path where 
<fr( + \z) is maximum at the saddle point. It can be easily found that the function has two 

real saddle points where ^ + ^(z) as well as its second derivative are real, but only one of them 
corresponds to a maximum of when the real axis is crossed perpendicularly. The latter 

saddle point is x c = exp [cosh" 1 , namely by using cosh" 1 (x) — hx[x + \/ x 2 — 1], 



x c = x + (j) + \Jx\(j) - 1 



(6.31) 



with 



x +(j) = - 
a 



f + bf + - 



(6.32) 



The constant phase contour T* which crosses the real axis perpendicularly at x c can be looked for 
in the form z*(6) — e x *^ +1 ° . It proves to be 



zJ6) = e cosh 1 (y*( M+ i( 



(6.33) 



where 8 g] — tt, 7r[ and 



y+(9; j) = u(8; j) cos9+ X /l + -u(0; j) + u 2 (6; j) 

a 



(6.34) 



with 



u(8;j) 



9 2 2f 



sin V a 

The contour V* crosses the negative real axis at the point z+(8 = ir) = —x+(j) with 



(6.35) 



x *(j) = - ex P 



cosh 



h 7T — 

a a 



(6.36) 



which lies on the cut ] — oo, — £>], because — x > — — exp [— cosh 1 (-)] . As a consequence, the 
unit circle can be deformed into the contour T* and a contour C{ _ , that goes around the cut 
] — oo, — x>] between the points — x*(j) and — x> in the clockwise sense (see FigHl 



dz 

|*|=i 27riz 



r $(+)( z ) _ 



dz 
— , — t 
27ri z 



r*<+'(z) , 



dz 
2tti z 



r$<+'(^) 



(6.37) 



On the other hand, in the integral involving (z) the unit circle can be deformed into a circle, 
minus the point on the negative real axis, with radius R that goes to infinity and a path around 
the cut ] — oo,— x>]. By using the parametrization z = Re 1 ®, with 6 ^ n, we get the following 

so that the contribution of 



large \z\ behavior: | exp [r$( '(• z )l I ~ ex P ~ T VaR cos(6*/2) 

|z|— >+oo L 

the integral along a circle of radius R vanishes in the limit where R goes to infinity. Consequently, 
(see Fig{5), 



| z |=i 



dz 
2T:i z 



r*<-'(*) _ 



]-«,-»>] 



dz 
27ri z 



r*'-'(z) 



(6.38) 



The crucial point is then to notice that the sum e T ® i+) ( z ) -). e r ^ { H z ) [ s an analytic function of z, 
which has no cut along the interval ] — oo, — sc>]. As a consequence, the integral along the contour 
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Figure 4: Deformed contour for e 1 "* + W which decomposes as a steepest-descent contour T* which 
goes through the saddle-point x c in the anti clockwise sense and a piece circumventing part of the 
cut ] — oo, — x>] in the clockwise sense. The other cut [x<, 0] lies inside the unit circle represented 
by a dotted line. 



*(+) 



with >{z) can be replaced by the opposite of the the same integral with & + \z) in 



place of )(z), and the equality (6.381 becomes 



dz 
=1 2tti z 



dz 
2ni z 



r* <+) (^ 



(6.39) 



where C- 



(-) 



is the contour which goes around the cut ] — oo, — x>] in the anti clockwise sense. 



i o 



Figure 5: Deformed contour for e T * ( W consisting of the circle at infinity, not represented here, 
and a path which circumvents the whole cut ] — oo, — sc>] in the clockwise sense. 



When the contributions (6.37) et (6.391 from the integrals involving respectively $(+) and 



are summed according to the definition (6.28 1 of K c (j; r), we get 



-) 



2K c (j;t) = ^ 



dz 
2ni z 



r* ( + ) (z) 



+ 



dz 
2i:i z 



r* <+) (z) 



(6.40) 



We stress that exp [r$^ + ^(z)] diverges when \z\ goes to infinity, except on the negative real axis, 
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so that the contour integral along the cut ] — oo, — ac*] cannot be closed at the point z = — oo. The 
expression (6.401 corresponds to integrate e T * t+> ( z )/(27ri z) along the contour in Figjfj] 




Figure 6: Contour of integration for e r$< W in the integral representation (6.40) of K c {j.t) 



On the contour C]_ 00i _ a ,^i , z = e A + 1CT7r where a = +1 if z is above the cut and a = — 1 otherwise. 



Since tj = n where n is an integer 



dz T<B (+) 



2ni z 



1)" / — e~ r3A sin(rVacoshA-6) . (6.41) 

Jinx, 7T V / 



The sign of this contribution changes for two consecutive values of j, but its absolute value is 
bounded, 



dz 
2ni z 



r$ ( + ) (z) 



1 

< — -e 

TTTJ 



— tj In x± 



(6.42) 



6.2.3 Large r behavior 

According to the saddle-point formula 



dz 
* 27ri z 



t— ^+oo 



X d 2 $( + )/dz 2 | z=2 



(6.43) 



By using the inequalities, — lnx* < — lnx c , derived from the expression (6.31 )-( 6.32) and (6.361 



and — j\nx c < $W(i c ), derived from (6.291, the bound exhibited in (6.421 implies that 



dz 
— , — i 
27ri z 



r$<+)(z) _ 



/ r $(+) 

o I e 



(x c ) 



(6.44) 



where o (e T ^^) denotes a function which decays faster than e T f^ when r goes to +oo. 



behaviors rt6.43 1 et (6.44) lead to 



Eventually, the definition (6.261 of K c (j;t) and the decomposition (6.401 together with the 



Ctj ( t )I T j integer T _^ c 



where 



and 



<t> i+) (x c ) = -j]n[x+(j) + y/x%G)-l] + ^b + ax + (j) 



1 1 



2V2^x c ^/ 'd?&+)/dz 2 \ z=Xa 



(6.45) 

(6.46) 
(6.47) 
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with x c = x + (j) + J x 2 + {j) — 1 and 



dHM/dz% =Xc 



(6.48) 



1 



The same argument can be performed for Ks An \j;r) defined in (6.271, with the result 

Sr 3 +A„| T3 intcgor (r) r ^ +oo A 7 x e ^ (+) (»e) with M = 
As a consequence, 

1 



and 



lim -ln[e- T p^Z T j(r)]=m 



lim. - In [e r p TJ s T j +A n (t)] = f(j) 



(6.49) 

(6.50) 
(6.51) 



with f(j) = —1 + j\sxp + $' + )(a; c ). By virtue of (4.56) Inp = In y j|jrf an d, according to the 
definitions in ([4T59]), cosh" 1 (^) = cosh" 1 ( /z|2 A B2 ) = In , / #§. Therefore /(j) reads 



/C?) = — 1 + j cosh" 



(V) 




-III In 





x+(j) + Jx\(j)-l 



+ \/b + ax + (j) 



(6.52) 



where x + (j) is defined in (6.321. 

Eventually, according to (4.481, P = n; r) is a finite linear combination of functions of n 
plus a finite increment An, which can be rewritten as 



V 9Au(tj + An;r), (6.53) 

Ara=0,l,-1 

and all functions <7a« prove to have the same "large deviation function" /(j) in the sense of 



definition (6.25 I 



.9a«(tj + Au;t) ~ ^a«(j, r)e 

r— >-+oo 
r J integer 



Therefore, by comparison with (6.23) and (6.241 we get 



fa 2 =f(j) and A(],t)=T ^ &An^An(j,T) 



(6.54) 



(6.55) 



An=0,l,-1 



where the expression (6.521 of f(j) indeed coincides with the result (6.101. 



6.3 Derivation by Laplace's method on a discrete sum 

As the reader may have noticed, the computation of the large deviation function via contour 
integrals is a bit tricky and clumsy due to the cuts, and relies on some compensations which are 
not totally obvious to foresee. 

In the case at hand it is possible to derive the large deviation function via Laplace's method 
applied to a discrete sum of non-negative contributions. We illustrate this briefly in the case of 



K. 



(An) 



37 



The key is an explicit formula for sin ^ry A + (z) j / y A + (z) as a Laurent series in z. From 
the symmetry z f> 1/z we can concentrate on positive powers of z. We start with 



sin(r v /A + (z) +oc 2k+1 

' rTT ^=E(^ T I)T[^)] fe (6-56) 
A + (^) fe=o v ; 



and expand [A + (z)] as a Laurent polynomial in z. 



2 V z / / ^ \1) l\m\(h-l-m)\ 



l,m>0 
l+m<k 



So for n > one gets (take / = m + n above) 



dz 1 s H-V A +(*)y _ _ ^ ,a^ + ". fc _ 2m _ n fc! 



E E f2*!+ 1)1(2) 



b 



,„,_! 2m /t , ; ^ ^ (2fc + l)!V2/ m!(m + n)!(fc-2m-n)!' 

(6.58) 

As t, a and & are > 0, this is a (double) sum of positive terms, and we are interested in the limit 

t ^ +00 n = Tj + An with Art e {0, — 1, 1}. (6.59) 

It is straightforward to check that in this limit the maximal term in the (double) sum is in the bulk 
(i.e. not for m — or fc = 2m + n) and such that fc and m scale linearly with r. One can use the 
Stirling approximation for all factorials and obtain the large deviation function straightforwardly, 
the most painful part of the computation being the location of the maximal term. We omit all 
details. 

7 Dependence upon typical time scales of the thermostats 

In this section we are interested in the influence of the typical time scales in the limit where the 
heat exchanges with one of the two thermostats become infinitely fast. We only explicitly consider 
the limit when 1^2/^1 goes to infinity, namely when the heat bath 2 exchanges heat with the two- 
spin system far faster than the heat bath 1 does. Indeed, the opposite limit when ^2/^1 vanishes 
involves similar calculations (where the roles of Qf and Q 2 are interchanged), and the results in 
both limits are essentially the same when they are stated in terms of quantities pertaining either 
to the "fast" heat bath or to the "slow" heat bath. The main results are summarized with the 
latter terminology in subsection |1.2| 



7.1 Stationary mean values in the infinite v^jvx limit 



The evolution of the probability distribution P(a\,a 2 ;t) is given in (4.401 where the time scales 



involve the parameter a = \/T — 4 A with A = V\Vi (1 — 7172). When the ratio of inverse time 
scales v 2 /v\ goes to +00, lim i , 2 / I/1 _ J . +00 \ (v x +v 2 ) [1 + a] = v 2 and lim I/2 / Wl _ J . +00 §(^1+^2) [1 - a] = 
1/1(1—7172), so that the probability distribution P(a±, 02] t) goes exponentially fast to its stationary 
value over the time scale l/[vi(l — 7172)]- Moreover, since 7* = Viji + V2J2, hhv 2 / ! , 1 _ i . +00 7* 



72, and, according to (3.6 1, the stationary probability distribution P s t(ci, 02) coincides with the 



canonical distribution at the inverse temperature P2 of the fast thermostat, namely 

lim Pst (01,02) = 7(1 + 720-1 £72)- (7-1) 

V2/V1— »+oo 4 
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As a consequence the heat capacities (Tj, T 2 )'s corresponding to a variation of the tempera- 
ture T a of bath a (for a = 1,2) and given in ( |3.22[ ), become \vca. 1/2 / v ^^ 00 (7^(71, T 2 ) = and 

lim y2 / yi _j. +00 Cg t (Ti, r 2 ) = C cq (T2) respectively, while the heat capacity C s t(7i,T2) correspond- 
ing to equal variations of both bath temperatures, and given at the end of (3.2.31, becomes 

nIr V 2 /fi-s-+co C st (Ti,T 2 ) = C C q(T 2 ). 

On the other hand the time scale of the mean currents of exchanged quantities is that of the 
slow thermostat. Indeed, according to (3.271, 



lim (j 2 )st = (71 - 72) 

I/2/l/l->+0O 



v\ Ae 



(7.2) 



and the thermal conductivity is determined by the typical time scale l/v\ of the slow thermostat. 
Similarly the housekeeping entropy flow ( |3.47| ) , which is equal to the opposite of the mean exchange 
entropy flow in the stationary state d^^S / dt\ si , becomes 



Ae 



lim cr hk [P st ] = Vi (71 - 72) {Pi - P2) ,, 
V2/V1— >+oo Z 



(7.3) 



7.2 Various probabilities in the infinite v<ijv\ limit 

The stationary probability that the heat received from the bath 1 be equal to Q\ = — niAe is 
given by (4.471 where 70 is to be replaced by lim^/^^+oo 7* = 72, with the result 



lim P s t(ni;t) 



lim e T [c ni (t) + s ni (r)] 
vz/vi— y+00 



(7.4) 



According to the definitions (4.351 and (4.36), the integral representation of c„(t) and s„(r) in 
the complex plane involve the discriminant A + (z). In order to discuss the dependence upon 
the inverse time scales V\ and i>2, it is convenient to rewrite the expression ( |4.22[ ) of A + (z) in 
terms of the parameters p + and p- defined in (6.151. The correspondence with A and B reads 
A = viv 2 (p+ + P-) and B = V^Vi (p + —P-), and 



- P- ) and B = v\V2 (p+ — p~ 
A+(z) = l + 2i/ 1 i7 2 



-(P++P-) +p + z+p- 



(7.5) 



Therefore lim,, 



1 v 2 /v 1 ^+oo A+(z) = 1 and the leading order of \{v\ + 1*2) [—1 + \J A+(z)] is merely 
\v\ [— (p + +P-) +p+z + p-(l/ z)]. Consequently e~ T c ni (r) and e~ r s„ 1 (r) have the same asymp- 
totic behavior and 



lim e 

u 2 /ui— S- + 00 



Kir) 



J\z\ = l 



dz 1 
, 2tti z n+1 



Grw{z; vit), 



where 



G RW {z]Vit) = exp<j ~^(P+ +p-)vit + 



1 


r 


Vit\ 








2 


z 





(7.6) 
(7.7) 



This expression can be interpreted as G R i V {z\vit) — Y^ti=-oo z "-FW( n ; namely G RW (z; U\t) 
is the generating function of the probability P RW (n;vit) for the continuous-time random walk, 
also referred to as the "randomized" random walk (see for instance page 59 of Ref.j^S]), which is 
determined by the Markov evolution equation 



dP mv (n; vit) _ 17 
dt ~ 2 



(p + + p-) P RW (n; v\t) +p + P RW {n - 1; v x t) +p-P RW (n + 1; v x t)] , (7.: 



and the initial condition P RW (n;t = 0) = <5„.o- By virtue of the identity which defines the 
generating function of modified Bessel functions I n (x), 



,i[p + z+p_i]yit 



+00 

E 



In (vity/p+p-) 



(7.9) 
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where 



2jt 



<m 
2V 



(7.10) 



Therefore the series representation of the expression (7.7) for the generating function G RW (z; vit) 
yields 



(7.11) 



Eventually the limit in (7.4 1 reads 



lim P st (ni;£) = P RW {ni]V\t), 

V2/V1 — >+oo 



(7.12) 

where the thermodynamic parameters of the thermal baths appear through the combinations 
of p+ and p- explicitly given in (6.181. On the other hand the probability that at time t 
the system has received a heat amount Q 2 — n 2 Ae from the bath 2 is given by (4.48). Since 
lim !/2 / !/1 _ ! . +0O e _T c„(r) = iP HW (n;i/ii) and lim„ 2/l , H+00 e~ T s„(r) = §P w (n; nt), 



lim P st (n 2 ;t) = 1 + l2 P RW {n 2 ; v\t) + 72 [Pw(^2 + 1; fit) + P R w{n 2 -l;vit)]. 

(7.13) 

We notice that at equilibrium, namely in the stationary state where j3\ — j3 2 — ft, p+ = p~ 



according to (6.181, and the limit (7.12) reads 



lim P eq (ni;i) = 7| ni | -(1 -7 )v x t 



Similar formulas hold for the three contributions in lim 



3 -i(l-7>lt 



(7.14) 



^2/^1— 5- + OC' 



P eq (n 2 ; t), which is derived from 



{7J3I. 

7.3 Interpretation: mean-field regime 

When the initial state is distributed according to the stationary measure, the probability that at 
time t the system is in a configuration where <j\a 2 is equal to ±1 and that the system has received 
a heat amount Qi — ~n 1 /S.e from the bath 1 is given by (4.451 and (4.46) where 7* is to be 
replaced by lim V2 / Vl ^. +0o 7* = 72, namely 



lim P st (a-102 = ±l,m;t) = - [1 ±72] lim e T [c ni (r) + s ni (r)] 

vi I vi— y+oo Z V2/V1— >+oo 



Comparison with (7.1 1 and (7.12) shows that the latter equation can be interpreted as 
P s t (oia-2 = ±l,ni;t) ~ P st (<J\(J 2 = ±1) x P st (n 1 ;t). 

V2/V1— >+oo 



(7.15) 



(7.16) 



This is a mean-field property: between two flips of the spin a±, the spin a 2 is flipped so many 
times by the thermostat 2 that, when the spin a\ is flipped again, the sign of (J\(J 2 is no longer 
correlated to its value when the previous flip of o\ occurred. Therefore the variation of n\, which 
is generated by the flip of o\ and the value of which is determined by the sign of a±a 2 , is no longer 
correlated to the sign which cti<72 had when the previous variation of ri\ occurred: the probability 
distributions of <j\ o 2 and n\ are independent from each other. 

On the other hand, the probability that at time t the system is in a configuration where ct\o 2 
is equal to 1 and that the system has received a heat amount Q 2 = n 2 Ae from the bath 2 has an 
expression given by the remark after (4.451 and (4.46). We get 



(I + 72) 2 

lim P st (a 1 a 2 = +l,n 2 ;t) = P RW {n 2 ;vit) 

v 2 /v 1 —y+oo 4 



P RW {n 2 + l;i^t), (7.17) 
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while 



lim P st (<ti<j 2 = -l,n 2 ;t) = — 72 - P R w(n 2 ;vit) + 72 P RW (n 2 - l;vit). (7.18) 

^2/^1^+00 4 4 



Comparison with (7.131 shows that there is no factorization similar to (7.161. In other words the 
variables n 2 and &\<j 2 are still correlated. 



7.4 Symmetry property specific to the probability of Q 1 

As can be derived from the integral representation ( |7.10| , the modified Bessel function I n (x) is 
an even function of n, J_ n = /„. Therefore the probability P RW (n;Vit) for the continuous-time 
random walk obeys the symmetry 



P RW (n;uit) p+ 
In — — - = n In — at any time t. 

P R w{-n;vit) p- 



(7.19) 



As a consequence, according to (7.12) and the relations (6.181, the probability distribution for the 
heat amount dissipated towards the slow bath, Qf = — niAe, obeys the finite-time symmetry 



In 



Pst(Qtt) 



P st (-Q?;t) »*A 



(7.20) 



However, by virtue of (7.13), there is no similar finite-time symmetry property for P s t(Q 2 ',t)- 



7.5 Cumulants per unit time for Qf and Q 2 



At any time, according to (7.7l, the characteristic function G RW {e \v\t) for the continuous-time 
random walk takes the very simple form 

G RW (e x ; vit) = exp[ta RW (X; v x )\ 

where 



(7.21) 



a RW (X; vx) 



- (p+ +P-) + P+e +P~e 



-x 



(7.22) 



On the other hand, according to ( 7.12 ), Yln 1 =—oo e ^ ni li m i/ 2 /^i-s.+oo 

P st (ni;t) = G RW (e x ;i/ 1 t). 

Therefore, the cumulants for Qf are given at any finite time by the formulae Knl = d q In G RW (e x ; v\t) / 'd\ q \\=o 
and the cumulants per unit time read for p > 



1 



J2p+1] 
t 711 



-w (P+ -P-) = ^(71 -72) 

1/2/^1—^+00 z z 



and for p > 1 



1 



= v (p+ = ^C 1 -7172)- 

£ 1*2 /vi -4-+oo Z Z 



-Pp] 



(7.23) 



(7.24) 



In the case of Q 2 , by virtue of (7.131 

+00 



]T e Xn *P st (n 2 ;t) 



1 



2/vi— >+oo 2 



[1 + 7 | + (1 - 72 2 ) cosh A] G w (e A ; i^t). (7.25) 



The cumulants per unit time of Q 2 coincide with the cumulants per unit time of Qf only in the 
long-time limit, namely for p > 



lim -/ 
t 



ir 11 = y(7r-7 2 ) 



(7.26) 
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and for p > 1 



lim 



1 



.[2p] 



^1 



t— >+oo i ™ 2 i, 2 / Vl — s.+c 



(1 - 7x72)- 



(7.27) 



When the system is at equilibrium 71 = 72 and the long-time behavior of all odd cumulants per 
unit time vanish, as already noticed in subsection |5.2| 

We notice that the previous results can also be retrieved directly from the expression for the 
generating function of the long-time cumulants per unit time. (The limits v 2 jv\ ~~ * +°° and 

t — > +00 do commute with each other.) The expression of a2(A) is given in (5.14), and, according 
to the relations after (6.151, it reads 



a 2 (A) = - I -(fi + v 2 ) + \j {v\ + v 2 ) 2 + 2^i v 2 



~(p+ +P-) + P+e x +p~e x 



(7.28) 



012(A) is a symmetric function of v\ and v 2 . In the limit v 2 jv\ — > +00 the generating function of 
the long-time cumulants per unit time becomes 



lim a%{\) = a RW {\;vx), 



(7.29) 



where u RW (\;v\) is given in (7.221. 

7.6 Long-time current distribution in the infinite v^jvi limit 



According to the definition (7.101 of the modified Bessel function 

/■27T Jf\ 

In=t j (cd)=J ^e**™ with g(6;j) = -ijd + a< 



(7.30) 



where a denotes some parameter. In the complex plane where the affix reads z = + iO' , g(z;j), 
the analytic continuation of g(0;j) is a periodic function of z with period 2ir when tj is equal 
to an integer. Therefore I n =tj can be rewritten as I n =tj = L- n n ~\(dz/27r) exp[tg(z;j)] and, by 
applying a saddle-point method to the latter integral, with a deformation of the initial contour 
in order to exhibit the constant phase path which goes through the saddle-point in the direction 



where it is indeed a maximum, as done in subsection 6.2 one obtains (with the relevant saddle 
point z c = -i ln[(|j|/a) + ^(j/a) 2 + 1]) that 



in=tj 



{at) 



t— >+oo 



2™V(j/a-) 2 



: exp 




+ 1 — |i| In 




The long-time behavior of a current density U(j; t) is given in terms of P(n — tj; t) by ( 6.22 !-( 6.23 1 
From the expression (7.11 1 for P RW (n = tj; Uit) with a = 17 ^Jp + p_ , we get 



(7.31) 



PZ,(j,t-,v 1 ) = A RW (j,t-,v 1 )e t f™^) 
with A RW (j,t; vi) = Xj^l-Kv^ty/p-p- + (j'/^i) 2 and 



(7.32) 




/rw0';^x) = - (P+ +P- )+j In 4 /— +| j I In y/p+P- +vi \ / ^ +P+-P-Hil m 
2 VP- V ^x 

(7.33) 

f RW {j] v i) is the large deviation function for the randomized random walk described by the Markov 
equation (7.8 1. We notice that, as predicted by large deviation theory and in particular the 



Gartner-Ellis theorem, f RW (j; v\) can also be retrieved as the inverse Legendre transform of the 
generating function of the long-time cumulants per unit time a RW {\;v\) written in (7.221. 
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The probability density of the cumulated heat current Qf jt in the long-time limit is given by 
(16.22b, (17.121) and 17.32b, with the result 



+00 



nas ( = x 1 J , = E 5 ( J - 7 ) M ™ y. *; 



n— — oo 



,i/swt)>i) 



(7.34) 



Similarly, according to (7.131, by an argument similar to that leading to (6.551, the probability 
ite< 

Q 2 



density of the cumulated heat current Q/xjt in the long-time limit is shown to read 



n 3 * \^r=r,t\ = 



J2 5 (j ~ -)tA 2 {j 1 t;v 1 )e 



,tf HW (j;ui) 



(7.35) 



with A 2 (j,t;i' 1 ) 7^ A RW (j, t; Vi). We retrieve that the cumulated currents — riiAe and Q 2 = 
n 2 Ae have the same large deviation function, and more precisely, 



lim fa 2 (J) = f RW I 

va/vi-t+oa \ Ae 



(7.36) 



We notice that the expression (7.331 for the large deviation function f RW (j; v\) agrees with the 
limit of the expression (6.16) for /(j), when v 2 jvi goes to infinity and j/v\ is fixed. Indeed, in 



the expression (6.161, which is valid when T\ ^ 0, the function Z(j) given in (6.17) is such that 

lim Z(J) = 4/^2 +P+P-> ( 7 - 37 ) 

j / v\ fixed 



while (v x + v 2 ) -1 + y/1 - 2V X V 2 \ V + +p-- 2Z(j)} ~ -v x ]p+ +P-- 2Z(j) 



8 Case where T\ = 0: pure energy dissipation towards the 
thermal bath 1 

8.1 Microscopic irreversibility 

When the temperature T\ of the colder bath vanishes, in the sense that /3iAe goes to infinity, the 



microscopic reversibility (2.2) is broken 



(— a, o"|W|er, a) = whereas (a, ojW| — a, a) 7^ 



(8.1) 



by virtue of the expression (2.9 1 for the transition rates when 71 = 1. In other words the thermal 
bath at zero temperature cannot provide energy to the system, i.e., it cannot flip the spin <j\ if 
the flip corresponds to an increase of the two-spin system energy. There is only energy dissipation 
towards the zero-temperature bath. 



When fa = +00, the ratio in the modified detailed balance (2.3 1 vanishes or is infinite when 



two configurations differs from each other by the sign of ci, and when the spin a± is flipped by 
the thermostat 1 the exchange entropy variation SSi H (C <— C) = —f3i8qi{C <— C) is infinite. All 
direct consequences of the modified balance are no longer valid. 



However the Markov matrix (2.10) of the configurations evolution is still irreducible (see the 



definition after (2.10 1 ) because histories such as 



(a, a) -> (cr, -a) -> (-a, -a) -> (-a, a) 



■2) 



correspond to a succession of flips with non-vanishing transition rates. Therefore, according to 
the Perron-Frobenius theorem, there still exists a single stationary distribution and in the latter 
distribution every configuration has a non- vanishing weight. The stationary probability given by 
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(3.6) is still a canonical distribution with an effective inverse temperature = (2/Ae) tanh 1 7", 
with 7° = V\ + V 2 j 2l 



lim P s t(<T 1 ,a 2 ) = -r[l + 7*crio- 2 ]- 

3iAe— 4 



According to (3.271, the mean current is finite, 



lim (j 2 ) st 

1 Ae— J'+OO 



,., N Ae 



(8.3) 



(8.4) 



Since j5\ = +00 the stationary exchange entropy flow (3.34) is infinitely negative in the stationary 
state 



lim 

lAe— >+oo dt 



(8.5) 



while the rate of entropy production, which has the opposite value in the stationary state, is 
infinitely positive, lim^Ae-^+oo dintS' SG /<^|st = +00. We also notice that the heat capacity with 



respect to a variation of the temperature 7\ from the zero value, C£ ] (Ti = Q,T 2 ), defined in p~19| 



vanishes according to (3.221 



As a consequence of the fact that the thermal bath at zero temperature cannot give energy to 
the system: Q\ = — niAe is necessarily negative and 



P(Qi,Q2\t)=0 if Qi > 0. 



.6) 



This can be checked on the explicit expressions of subsection 4.2 as follows. The probability that 
the system is in configuration (a%, &%) at time to = 0, in configuration (a[, a' 2 ) at time t and receives 
Qi = — nAe and Q 2 = (n + An)Ae during the time interval [0,t] is cr 2 |U(rt, n + An; t)\ai, <j 2 ). 
According to ( 4.37 1-( 4.38 1 where 71 is to be set equal to 1, the latter matrix elements involve 



the functions c n (r) and s n (r) defined in (4.35) and (4.36). When T\ vanishes, 71 = 1, A — B = 
v\V 2 [\ — 71) (1 + 72) also goes to zero and, by virtue of (4.22 1, lim^Ae-y+oo A + (z) = 1 — 2A + 2Az. 
As a consequence, cosh(r-\/ A + (z)) = ^ p=0 [l/(2p)!]r 2p (l — j4 + ylz) p contains only positive powers 
of z, and so does sinh(r-\/ A + {z))/ A + (z) . Consequently c„(r) and s„(r) vanish for n < and 
lim ) g 1 Ae-i-+oo(o'i,cr2l l L | ( n 7" + ^";*)|o'i,a'2) = for any n < 0. 



8.2 Long-time behavior 

The explicit values of the infinite-time limit for the cumulants per unit time of the heat and Q 2 

When Xi vanishes, 71 = 1, A tends to A a = V\V 2 {\ — 72), 



5.2 



can be calculated as in subsection 
and A — B vanishes. Then the expression (5.141 of a 2 (X) is reduced to 



lim a 2 (X) 

Pi Ae— >+oo 



(vi + v 2 ) + 4/ (vi + u 2 ) 2 + 2viv 2 {l - 72) 



■7) 



The expressions of the first three cumulants can be retrieved by setting 71 = 1 in the expressions 



(5.16) 



As in section [6] the large deviation function can be derived either as the L egen dre transform 

for P Po (n 2 ;t) 



6.2 



of a 2 (\) or by a saddle-point method similar to that performed in subsection 
given by ( 4.48[ ), which also provides the amplitude of the probability, or it can be retrieved di- 
rectly by taking the limit 71 -> 1 in the expressions ( 6.13 )-( 6. 14 1 for the large deviation func- 
fg 2 (j), as follows. In the limit where T\ vanishes, so does A — B, and Y(j) tends to 



tion 



Y°(j) =3 3+ \/f + 1-2,4° where A = v x v 2 {\ - 72), while, according to (|3.27[) and K3§ 



(j)st = ViV 2 {l - 72)- 
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As a result, 



lim Jq 2 (J) = -oo 



(8.9) 



lim f Q2 (J) = 



^1 + ^2 



(-1 + [1 + ln0)° t ] J + ^ + 1-20)0 - j In J [j + ^ + 1-20)°] ) 



where J = [2/(i/i + i/ 2 )] J/Ae 



8.3 Limit where v^jvx becomes infinite 
8.3.1 Finite-time behaviors 

The discussion can be performed along the same lines as in section [7] We only point out the 
features which are qualitatively different when 71 = 1. The discrepancies are due to the fact that, 
according to (6.151, when 71 = 1, though p + remains finite, 

p_ = 0, (8.10) 

and the random walk process associated to the variation of the heat amounts Qf or Q2 can have 
only positive increments. As shown below, a Poisson process shows off as a randomized random 
walk which can proceed only in the sense of increasing positive 712. 

Since p_ vanishes, while p + = (1— 72), the expression (7.5) for A + (z) becomes lim ( 3 1 Ae-s.+co A+(z) 
1 + 2l'iF2(l — 72) [ — 1 + z] and the generating function G RW (z; u\t) which appears in subsection 

[—l+z]p + v x t 



7.2 (see (7.7)) is to be replaced by 



G 



Pois 



>; v\t) 



with 



^i = ^(l-72>i= lim (j)° t = {j)° st * 

•v + OO 



(8.11) 



.12) 



Gp i S (z; v\t) is the generating function Gp i S (z; v\t) = J2 n =o z ™ft'ois(«; v\t) for the Poisson process 
ruled by the Markov evolution equation 



dP Pois (n; v\t) 
dt 



= p+vi [-Ppoi S (n; v\t) + Pp \s(n - 1; Uxt)} , 



and the initial condition Pp i S (^; t — 0) = <5 ri! o- The solution reads 

"Poisl n j v v-) — ■ e 



(8.13) 



(8.14) 



Eventually the probability that at time t the system has dissipated a heat amount Qf — niAe 
towards the bath 1 at zero temperature reads 



lim lim P st (m;t) = Pp i s (ni; V\t). 



(8.15) 



A formula similar to (7.131, where P HW (n 2 ; vit) is to be replaced by -Pp i S (™2; ^lt), holds for the 
probability that at time t the system has received a heat amount Q2 = 77,2 Ae from the bath 2. 

8.3.2 Long-time behavior 

The long-time cumulants per unit time for the heat Q 2 can be obtained by noticing that, when 



v\jv2 vanishes, the expression (8.7) for lim^A 



e— f+oo 



012(A) becomes 



lim lim o 2 (A) = a Pois (A; v x ) , 

Vafvi-t+aa /3iAe->+oo 



(8.16) 
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where 



(A;^) = ^(l 



12) 



-1 



is the cumulant generating function for a Poisson process with average v\p + = ^(1—72). 

■everse order: by virtue of ( 7.29 ) He 
apois(A; v\) . The long-time cumulants per unit time are all equal 



(8.17) 
We notice 



that the two limits can be taken in the reverse order: by virtue of ( 7.29 ) lim^ a 
(A; it 



+00 lim^/^^+oo a 2 



lim, 



lAe- 



lim - 

t— >+oo t 



lim 



lim 



V2/V1— ¥+00 /3iAe— >+oo 



(1-72). 



(8.18) 



The large deviation function can be retrieved 

- either as the Legendre transform of ct 2 (\) = ctp i s (A; v\ 

- or by a saddle-point method applied to the expression 
that performed for P RW (n; fit), with the result, 



14) of Pp i s (n; fit) and similar to 



£pois(n; vit)\ t 



1 



with 



\tj integer ^ +QO y^rfj 

/pois(i;^i) 



ot/pois(j';^i) 



<j>r 



(8.19) 



.20) 
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where (j)g t * is defined in 

- or directly by taking the limit v 2 jv\ — > +00 with jjv\ fixed in the expression for the large 
deviation function lim^Ae-^+oo /q 2 (^) given in (8.9) (and by noticing, that (j)!? t is of order ^ 

while ~f is of order (^-) 2 x (^) )i 

- or by taking first the expression for J > of lim l , 2 / I , 1 _ i . +00 fQ 2 { J) given by (7.361 and (7.331 
and then taking the limit /3iAe — » +00, namely p_ — > 0. 

Eventually, the large deviation function takes the simple form 



lim lim fa 2 (J) = 

i/ 2 /vi->+ooftAe-++oo j<0 



-OO 



.21) 



lim lim f Q2 (J) = f Pois (J) 

V2/V1— ¥+00 piAe— i-+oo ]>0 



1 

Ae 



where (J2 



,0* 



(z/i/2)(l- 72 )Ae. 



The expression of {S 2 ) s t c & n be interpreted as follows. Since the bath 1 is at zero temperature, 
the spin crj may be flipped only when it is opposite to the spin a 2 . Moreover, since i/ 2 > vi, once 
the spin 1 has been flipped so that <j\a 2 = 1, on average the spin 2 is flipped a great odd number 
of times with a net energy transfer Ae from the heat bath 2 until the spin o\ is again flipped with 
an energy transfer Ae to the heat bath 1 so that u\g 2 — 1 again. As a consequence the mean 
current energy is equal to Ae times the typical inverse time v\ between two possible flips induced 
by the thermal bath 1 times the probability that a 2 is opposite to o\. namely (1 — 7 2 )/2. 



8.4 Limit where v<ijv\ tends to zero 

In the reverse limit where v 2 <C v\ , the roles of the two heat baths in the discussion of section [7] 
are interchanged. The slow thermostat is the heat bath 2 and the evolution of Q 2 is a Poisson 
process with the kinetic parameter v 2 . 

The stationary mean heat current received by the system is now {Sf 2 ) s t — (^ 2 /2)(l — 7 2 )Ae. 
The interpretation of the latter expression is the following. Since v 2 <C v\, as soon as the spin 
a 2 is flipped to a value opposite to the spin o\ with an energy transfer Ae from the heat bath 2, 
the bath 1 flips the spin a\ so that <j\o 2 = 1 and an energy Ae is transferred to the heat bath 
1. The next flip can be only a flip of the spin a 2 and its probability per time unit is the value of 
the transition rate of spin 2 when a\u 2 — 1, namely (l/2)z^ 2 (l — 72). As a consequence the mean 
current energy is equal to Ae times (v 2 /2){l — 72). 
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9 Thermal cycles 



As recalled in the introduction, part of the physical relevance of the two-spin system is as an 
idealized mesoscopic thermal machine, with heat flowing from the high temperature reservoir to 
the low temperature reservoir in average. This flow of heat results from thermal cycles made by 
the system. After each thermal cycle, the two spins have returned to their original state, but an 
amount of heat 2Ae has been transferred from heat bath 2 (the hot bath) to heat bath 1 (the cold 
bath). Thermal fluctuations do occur however, and with this interpretation two questions come 
naturally. What is the distribution of the time needed to make a thermal cycle ? What is the 
probability that the machine will perform a thermal cycle in the wrong direction ? 
The graph showing the possible transitions in the two-spin system looks as follows 

(+,+) (+,-) 

t t • 

(-+) (-,-) 

This graph looks like a square, i.e. a cycle with 4 edges. Most of what we shall have to say 
applies equally well to a Markov process with a finite number of states and whose associated 
transition graph is a general cycle. As the probabilistic reasoning is more transparent in this more 
general framework, we shall devote a separate section to it for completeness. It is likely that the 
forthcoming analysis has already been performed (more than once) in the literature, but we have 
not found it. 



9.1 Markov processes with a cyclic transition graph 

We label the n > 3 configurations as 1, ■ ■ ■ , n, and identify configuration m with configuration 
to + n in all subsequent formulae. We view the process as the motion of a particle along the 
cycle, jumping from time to time from a site to one of its two neighbors. We choose arbitrarily 
an orientation of the cycle. A jump from to to to + 1 (resp. to — 1) is said to be clockwise (resp. 
anti clockwise). We let c m be the transition rate from configuration to to configuration to + 1 and 
a m be the transition rate from configuration to to configuration m — 1. All other transition rates 
vanish. We set b rn = a m + c, m . If the particle sits at to at time t, the probability that the next 
jump will be clockwise (resp. anti clockwise) is c m /b m (resp. a m /b m ). By saying that the graph 
associated to the Markov process is a cycle we mean that all a m 's and c m 's are > 0. To be totally 
explicit, with the conventions of this article, the generator of the Markov process looks like 
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C n -3 


-b n -2 


On-1 















Cn-2 


-b n -l 


a n 
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Cn-1 


-b n 



Starting at some arbitrary reference configuration, each later visit to this configuration defines 
an integer, namely the algebraic number of times the cycle has been traversed, i.e. the homotopy 
class of the trajectory between the two passages at the reference configuration, or the winding 
number. The variation of winding number between two successive visits to the reference configu- 
ration belongs to {—1, 0, 1}. We can refine the definition of winding to have it defined at all times, 
i.e. as a process W t , via the following trick: nW t is just the number of clockwise jumps minus the 
number of anti clockwise jumps that have occurred up to time t (included). Note that Wt is an 
integer if and only if the positions of the particle on the cycle are the same at time and t, and 
then Wt is simply the previously defined winding number. 
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9.1.1 Forward and backward thermal cycles 

The easiest question to answer is whether or not the winding number will ever reach ±1 ? The 
(strong) Markov property is the crucial ingredient. 

Let II~ be the probability that the winding number of a trajectory started at m ever reaches 
the value — 1/n. Then by the Markov property 

n m = T i~ t n m+1 n m . (9-2) 

Dm "m 

The meaning of this equation is clear : either the first jump is anti clockwise (probability a m /b rn ) 
and the winding number reaches its target —1/n or the first jump is clockwise (probability c m /6 m ), 
and then the particle has "lost" a winding 1 /n so it has to go from m + ltom with winding number 
— 1/n to compensate (probability II~ +1 ), and take a new chance. 

Luckily, we do not need to solve the full system: the probability that starting from m the 
winding number reaches —1 is, by the Markov property again, UlLn 1 ^m-i> which is independent 
of m. So we denote this probability simply by II = Jl/Lo 1 ^Kn— V We rewr ite the previous equation 
as 

a m (l - U m ) = c m n m (l - n m+1 ). (9.3) 

A first consequence is that if II^ = 1 for some m then also II~ +1 = 1 and so on, so that all II~ 's 
are equal to one, and the probability to reach winding number —1 is unity. On the other hand, in 
the case when no LI~ equals one, 



n = ft — (9-4) 



Indeed, in (9.3 1 take the product over all m's in the cycle and simplify both sides by rim=i(l — 
n m) ^ to get n'Ui a m = Ilm=i c m Ilm=i n ~ > Le - n_ = U m =i a m /c m . As IT is a probability, 
this is possible only if []™ =1 a m /c m < 1. 

We have already proved the following : if Y[ m =i a m/ C m > 1 then II - = 1, i.e. winding 
number —1 is reached with probability one. We could reproduce the above argument with 11+ , 
the probability that the winding of a trajectory started at m ever reaches the value 1/n, and 
n + . We would get : if Y[ m =i a ™/ C m — 1 then II + = 1, i.e. winding number 1 is reached with 
probability one. In particular, if Om=i a m = Ilm=i c ™ then the probability to reach winding 
numbers —1 and 1 is unity, and by the Markov property, the probability to reach any winding 
number an infinite number of times is also unity: the winding number W t will oscillate and take 
arbitrarily large positive and negative values as t — > +oo. 

To deal with the case Hm=i a m ^ c mi w e need a deeper result which we shall not prove 

here: the ergodic theorem for a finite state Markov process, which, stated informally and adapted 
to the case at hand, says that Wt/t will be close to (Wt) s t/t with probability close to 1 at large 
times. We shall "prove" below a formula for (W t ) st : 

(w^^ niUcr.-niUa^ (95) 

P n -i[a.,c) 



where P n _i(a.,c.) is the principal minor of the generator ( |9.1| of the Markov process (note that 
by construction the determinant of the generator is 0) . It is a homogeneous polynomial of degree 
n — 1 in the a m 's and c m 's given explicitly by 

n n 

p n _i(o.,c, ) = ^2 ^2 Yi a ™+i n Cm+k > ( 9 - 6 ) 

1 = 1 m=l l<j<l l<k<n 



where as usual an empty product stands for 1. We observe simply that the denominator in (9.5 1 is 
a sum of positive terms so that the sign of (Wt) s t/t is that of the numerator. The above discussion 
implies that 

n n 

n- = min{l, J] n+ = min{l, J[ ^}. (9.7) 

m=l ° m m=l a ™ 
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We observe that the quantity A = — lnmin{II~, II + } plays an important role in the context of 
chemical reactions, where, after multiplication by k B T, it is called the affinity of the cycle properly 
oriented. If A = — logII~ (resp. A = — logll" 1 ") the reaction has a tendency to evolve clockwise 
(resp. anticlockwise). 

Coming back to the physical interpretation, we may view the cycle as a mesoscopic thermal 
machine, and imagine that when the winding number changes by one, some reference heat bath 
has collected one unit of energy. This is what happens in the simple two-spin system, the "unit" 
being 2Ae. For arbitrary n, a microscopic implementation of this behavior is not so obvious, but 
is not needed for the discussion either. By the ergodic theorem for a finite state Markov process 
as recalled above, W t /t will be close to (W t ) s t/t with probability close to 1 at large times. So the 
behavior of the thermal machine is deterministic at large times. But fluctuations may occur, and 
e~ A is the probability that a time exists at which the net heat transfer is opposite to that of an 
expected working cycle, i.e. the thermal machine has performed the equivalent of a working cycle 
in the wrong direction. By the Markov property, e~ kA , k = 1, 2, • • • is the probability that a time 
exists at which the net heat transfer is —A: times that of an expected working cycle. 



9.1.2 Fluctuations in the time it takes to make a thermal cycle 

We can now come to the second question, namely what are the fluctuations of the time it takes 
to make one cycle in the direction the machine is built for? 

As usual, the object that satisfies simple equations is a Laplace transform. We assume that 
llm=i c m > Ilm=i a "i> i- e - tnat winding number grows in average. We denote by T m the random 
time it takes, starting from to to reach winding number 1/n. Our real interest is the random 
time T it takes to reach winding number 1, and by the (strong) Markov property, T is distributed 
as a sum of n independent random variables each of which is distributed like a T m . We write 
F m (X) = (e~ XTm ) m where (• • -) TO is expectation with respect to trajectories starting at to. This 
is a bit redundant here, but we want to stress that T is a cyclic invariant so that we may write 
F(X) = (e~ AT ) because the expectation with respect to any initial probability distribution gives 
the same result. 

The Markov property says that F — Ilm=i ^m- The F m 's satisfy 

F - Cm I " m F ,F (9 81 

m + A m + A 

The interpretation of this equation is analogous to that for II* . If the particle is at to at some 
time, it waits an exponential time with parameter b m and then jumps to to + 1 (resp. to — 1) with 
probability c m /b m (resp. a m /b m ). The computation of the Laplace transform of the waiting time 
gives 

/ + °°^ m e^* e - At = -^-, (9.9) 
Jo b m + X 

which multiplied by the jump probability c m /b m (resp. a m /b m ) yields the prefactors above. Then 
writing T m as a sum of the exponential waiting time to leave to and some other (independent of 
the waiting time at to by the Markov property) random time, we note that if the jump is to to + 1 
this random time is 0, whereas if the jump is to to — 1 this random time is, by the (strong) Markov 
property, distributed like the sum of two independent random times, one distributed like T m _i 
and the other like T m . 

A standard trick to deal with the quadratic equations for the F m 's is by a Riccatti transform, 
to linearize them. Then F(X) acquires an interpretation as an holonomy. So we write F m = 
Gm-i/Gm for to = 1, • • • , n and observe that F — Go/G n : whereas the sequence F m is periodic 
by construction, the sequence G m is not, and F is the holonomy along the cycle. With this in 
mind, we set G_i = FG n -\. Then we define 

F m = ( ( b -+^)/ a ™ -cm/a m ) ( 9 1Q ) 
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and check that (9.8 1 turns into the linear equation 



(g™-i ) — ( G m ') 



(9.11) 



valid for m = 1, 



Setting F = Fi ■ • • F w , one finds by iterating the above formula that 



F 



so that F is an eigenvalue of the transfer (or Bloch-Floquet or 
matrix F, i.e. a solution of 

F 2 - F Tr F + Det F = 0. 



(9.12) 

depending on the community) 
(9.13) 



This formula shows clearly that F is a cyclic invariant, because Tr F and Det F are. By continuity, 
F(0) = 1, and we shall see shortly that this allows to choose the right branch. 

Equation (9.121 allows to express F n — G n _i/G n in terms of F and the matrix elements of F. 
A moment thinking shows that the same formula expresses any other F m in terms of the same F 
and the matrix elements of the matrix obtained by applying a cyclic permutation of order m to 
the factors defining F. 

The determinant of F is easily seen by multiplicativity to be e A = Ilm=i c m/o,m (which by the 
way does not depend on A), so the affinity has also something to say on F. 

The trace of F can be computed in any specific case, but it is complicated even for our simple 
two-spin system. We conclude this short digression in the world of general cycles by computing 



F to first order in A. This will allow us to give a "proof" of formula (9.6) for (Wt) s t- We observe 
that F m (A) = F m (0) + X/a m (§g), so that 

TrF'(O) = Tr(f2 f i(°y-- f ™-i(°)^( J q V m+1 (0) • • • F„(0) ) (9.14) 



n 

— ' ' ' MO)Fi(O) • • • F m _i(0)] u 



(9.15) 



where [•••]n stands for the top left corner matrix element of the matrix product. Expanding 
d9T3| to first order in A we get 1 - Tr F(0) +Det F = and [1 - Det F] F'(0) = Tr F'(0). The first 



equation implies that the right branch for (9.131 is (under our assumption that A > 0) 



F 



Tr F - ^/(Tr F) 2 - 4Det F 



and the second then gives 



(T> = 



TrF'(O) 



Det F - 1 

By computing Tr F'(0) we shall now show that 



(T) = 



P n _i(a.,c.) 



n 

11/; 



n 

11/; 



(9.16) 



(9.17) 



(9.18) 



This is nothing but the inverse of the value announced in (9.5 1 for (Wt) s t/t- That the relation 
(T) (Wt) st /t = 1 should hold is intuitively clear: (Wt) s t A is the average of the growth of the winding 
number per unit time, and (T) is the average time it takes to increase the winding number by 
one unit. Intuition is not proof, but modulo that, we shall have given a proof of the formula for 
(Wt) s tA- J us t note that in (T), the expectation is with respect to any initial distribution, while 
the average of Wt is exactly proportional to t only if expectation is taken with respect to the 
stationary measure. 

We sketch the derivation that (JlILi fl i)Tr F'(0) = P„_i(a.,c), which gives immediately the 
announced formula for (T) . We first note that F m (0) has ( } ) as right eigenvector with eigenvalue 
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1 whatever the value of m, so that the same is true for any product of F m (0)'s. The generic 
matrix with this property can be parameterised by K(x, y) = ( x ^ Zy )• We observe that F m (0) = 
K(c m /a m , 0) and a simple recursive computation shows that K(xi, 0) • • • K(xk, 0) = IK^fc+a^iCfc-i + 
■ ■ ■ + XkXk-i ' • • X\,Xk + XkXk-i + ■ ■ ■ + XkXk-i ■ ■ ■ X2). In particular, the top left matrix element 
of F m+1 (0) ■ • • F„(0)Fi(0) ■ • • F m _i(0) is 

-. Cm — l Cm — 1 ' ' ' C\ C m —\ • • • CiC m _|_l C m _i ' ' ' CiC m _|_i ' ' • C n , . 

1 H r • • ■ 1 r • • • H ■ (y.iyj 

«m-l O-m-1 ' ' ' a l Qm-1 ' ' ' Ql a m+1 "m- 1 ' ' ' &1 a m+l ' ' ' a n 



Multiplying by (Jl"=i a i)/ a m and summing over m one recovers the formula (9.6) for _P n _i(a.,c. ) 



9.2 The case of the two-spin system 

We apply the formulae of the previous subsection to our concrete model. We choose the cyclic 
order 

(+,+) (+,") 

t I , (9.20) 

(-,+) <" (-.") 

with corresponding a; c transition coefhcients (remember c is for the clockwise transition and a for 
the anticlockwise one) 

!/i(l-7i)/2;^ 2 (l-72)/2 <-> i/ 2 (l + 72)/2;i/i(l+7i)/2 

t t . (9.21) 

! y 2 (l + 72 )/2;j/ 1 (l + 7l )/2 o ^(1 - 7l )/2; i/ 2 (l - 72 )/2 

The meaning of this diagram is that, starting at the upper-left corner (+, +) for instance, the rate 
from (+,+) to (+, — ) is the clockwise coefficient ^ 2 (1 — 72 )/2, while the rate from (+, +) to (— ,+) 
is the anticlockwise coefficient vi(l — 7 i)/2. 

The product Y[m=i a m/ C m is readily evaluated to be e~ A = e _2 ^ 1_ ^ 2 ^ Ae . Note that this 
quantity has a purely thermodynamic interpretation : this is due to the modified detailed balance. 
Under our standard assumption j3% > /3 2 , we see that in average the thermal machine works in 
the clockwise direction, and the probability that at some time the net heat transfer is — 2fcAe 

(fc = 0, 1, • • • ) iS e -2M/3i-ft)Ae_ 

Using the formula; of the previous subsection, a closed formula for (e _AT ) is not difficult to 
write down, but it is very complicated and not particularly illuminating. We content to give the 
first two cumulants: 

CO = (- + -) ^4r> ( 9 - 22 ) 

(T 2 > - (T) 2 



V\ V2j 7 1 - 7 2 

_ 1 1 8 

( 7 i - iiY 



(1-7172) ( 4 + +(2- 7l 2 - 7 |) ' 



1 V $J " " v\vi 



As expected, the cumulants of T diverge when 7 i — 72 \ : for 7 x = 72 there is no net heat 
current and the winding number is just a simple symmetric random walk. 

We conclude by observing that the formulas for the normalized random variable 



are slightly simpler: writing 



f = T( 7l - 72 ) 2 iW2 (9.24) 



log(e-^) ee g ^A fe , (9.25) 



2(71-72) m ' f^kl 

the pkS are polynomials in 7 i, 72 ,^i,^ 2 , with integral coefficients, homogeneous of degree k in 
^11^2 by dimensional analysis, non-homogeneous but of degree 2k — 2 in 7 i, 72 and symmetric 
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both under the exchange of v\, and under the exchange of 71, 72- It is easy to check that in the 
limit when "fi = 72 = 7 the formulas simplify dramatically. In fact, setting v = v\ + ^2, one has 

H» (9.26) 

7i.72->7 i(7l — 72j (.1 — 7 J 

a result reminiscent of the theory of the simple symmetric random walk. 
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